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PREFACE 


Many  practical  problems  in  industrial  and  social  planning  require  op¬ 
timal  decisions  to  be  made  periodically  through  time.  Linear  programs,  called 
dynamic  linear  programs,  can  often  be  formulated  to  model  the  requirements 
of  these  decision  processes.  These  programs  are  generally  quite  large  and 
difficult  to  solve.  The  search  for  efficient  methods  in  finding  their  optimal 
solutions  has  been  a  major  topic  in  operations  research  for  the  past  25  years. 

Often  the  problem  modeled  by  a  dynamic  linear  program  involves  un¬ 
certainties  which  can  complicate  and  exponentially  increase  the  program’s 
size.  There  may  be  several  possible  outcomes  in  the  future,  but  determinis¬ 
tic  linear  programs  usually  only  consider  the  most  likely  outcome.  In  this 
dissertation,  we  present  methods  for  solving  the  stochastic  dynamic  linear 
program,  the  dynamic  linear  program  with  uncertainties  explicitly  included. 

Our  methods  take  advantage  of  the  program  structure.  Dynamic  linear 
programs  are  characterized  by  a  staircase  structured  coefficient  matrix,  in 
which  non-zero  elements  appear  only  in  blocks  along  the  diagonal  or  adjacent 
to  the  diagonal.  This  structure  makes  many  efficient  techniques  possible.  We 
will  show  that  the  stochastic  model’s  specific  structure  can  lead  to  additional 
procedures,  and  that  these  procedures  may  improve  upon  complicated  “brute 
force”  solution  methods. 

We  begin  in  Chapter  I  by  presenting  sufficient  conditions  for  a  deter¬ 
ministic  problem’s  optimal  solution  to  solve  a  stochastic  problem.  The  second 
chapter  discusses  the  difficulties  involved  in  using  deterministic  solutions  in 
general.  We  also  explore  the  possibilities  for  combining  separate  deterministic 
solutions  and  give  examples  of  problems  that  require  the  stochastic  dynamic 
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linear  program  to  be  solved. 

Chapters  HI,  IV,  and  V  present  methods  for  solving  the  full  stochastic 
program.  The  first  method  follows  from  the  decompositions  approach  to 
large-  scale  programming.  The  next  two  methods  employ  different  large-scale 
structured  programming  techniques,  in  which,  the  basis  is  partitioned  but  not 
completely  decomposed. 

Chapter  VI  demonstrates  that  the  methods  we  present  are  actually 
dynamic  programming  approaches.  They  only  differ  in  their  strategies  for 
approximating  the  optimal  state  space  solution  at  each  stage. 

In  the  final  chapter,  we  present  some  computational  results  for  our 
algorithms  and  discuss  potential  areas  of  applications.  We  also  state  our 
conclusions  on  the  use  of  stochastic  dynamic  linear  programs  and  suggest 
areas  of  future  research. 

In  this  dissertation,  we  use  standard  mathematical  notation.  More  specific 
notational  conventions  are  defined  in  the  text.  Within  each  chapter,  we  refer 
to  equations  and  propositions  by  their  order  of  presentation  in  that  chapter 
(eg., equation  (12)).  In  reference  to  equations  in  other  chapters,  we  prefix  the 
equation  by  the  chapter’s  roman  numeral  (eg.,  (11.12)). 
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CHAPTER  I 


Deterministic  Solutions  for  Stochastic 
Dynamic  Linear  Programs 


1.  Introduction 

Dynamic  linear  program  models  have  been  formulated  for  many  different 
practical  situations.  When  these  models  involve  uncertain  quantities,  the 
solution  of  the  resulting  stochastic  dynamic  linear  program  can  be  very 
difficult.  Under  some  circumstances,  however,  an  associated  deterministic 
problem  can  be  stated  that  is  easier  to  solve  and  can  be  used  as  the  “solution” 
to  the  stochastic  program. 

In  this  chapter,  we  will  state  conditions  that  imply  that  this  deterministic 
solution  is  in  fact  optimal  for  the  stochastic  dynamic  linear  program.  We 
call  a  solution  “deterministic”  it  solves  a  program  that  does  not  allow  for 
any  uncertainty  in  the  program  parameters.  We  also  will  use  myopic  solution 
to  refer  to  a  solution  of  a  program  in  period  t  that  does  not  make  use  of 
information  from  periods  after  t.  A  solution  that  considers  uncertainty  in 
the  future,  is  called  stochastic. 

In  our  development  of  an  optimal  deterministic  solution,  we  first  present 
the  basic  multi-stage  model  and  the  various  approaches  taken  for  its  solution. 
Also,  in  Section  2,  we  introduce  some  terms  and  notation  that  appear  in 
the  following  chapters  and  discuss  the  value  of  having  information  about 
the  random  variables.  To  do  this,  we  present  inequalities  that  measure  the 
superiority  of  the  stochastic  solution  over  a  deterministic  solution. 

The  chapter  concludes  in  Section  3  with  the  description  of  conditions 


for  an  optimal  deterministic  solution.  We  also  give  an  example  of  a  prob¬ 
lem  in  which  a  deterministic  solution  is  optimal  and  illustrate  how  added 
complications  can  necessitate  a  stochastic  solution. 

3.  The  Multi-Stage  Problem 

The  dynamic,  or  multi-stage,  linear  program,  to  which  we  shall  refer, 
has  the  following  form: 


min  =  cjiEi  •  •  •  -f  er^r 

subject  to  A\x i  =  b i, 

=  fa  (DLP) 


—Bt—iZt—i+AtZt  —  fr, 

zt  >  0  for  all  t, 

where  zt  £  Sfn*  (n^-dimensional  Euclidean  space),  £>i  £  S’7*1,  ft  €  &"**,  and 
the  vectors,  ct,  and  matrices,  At  and  B%,  are  dimensioned  to  conform. 

For  a  DLP  problem,  the  right-hand  sides,  ftl  are  given.  For  a  situation 
in  which  ft  is  random,  the  DLP  becomes  one  possible  program  out  of  the 
possible  outcomes  for  ft.  In  our  analysis  here,  we  will  not  let  any  of  the 
other  quantities  be  random,  so  ct,  At,  Bt,  and  b i,  will  be  assumed  known. 
Moreover  it  is  assumed  ft  and  fa  are  independent  t 

We  can  also  view  DLP  as  an  optimal  control  problem  by  assuming  more 
structure  for  the  matrices.  In  this  case,  we  would  have  zt  =  ( Vt ,  u«),  where 
yt  is  a  state  variable  and  u*  is  a  control  variable,  and  we  would  partition  the 
matrices  and  vectors  as 


At  = 


2 


(1) 


X 


and 


Bi  = 


(2) 


This  formulation  includes  transitions  from  state  to  state  according  to 


— A-tVi  —  Btut  —  yt+ 1 
and  interperiod  requirements  of 

Gt+iXt+i  +  =  f«+i- 

The  random  vector  x  represents  unknown  state  to  state  transitions  and  fu¬ 
ture  requirements.  Although  some  computational  efficiency  may  be  afforded 
by  the  special  structure  of  this  model,  we  will  restrict  our  discussion  to  the 
general  case  of  DLP. 

The  first  approach  we  consider  for  DLP  is  to  solve  it  for  all  possible 
ft  and  then  take  the  expected  value  of  z\  as  the  measure  of  cost.  This 
approach  requires  perfect  information  of  the  outcomes  of  all  future  events 
and  is  known  in  the  literature  as  the  “wait-and-see”  solution  (see  Madansky 
[38]).  It  would  be  the  optimal  solution,  if  one  could  somehow  wait  and  see 
until  the  end  of  the  planning  horizon  and  could  beforehand  make  decisions 
based  on  what  will  occur.  Obviously  such  a  perfect  information  solution  is  not 
implementable.  When  averaged  over  every  possible  ft,  it  provides  a  measure 
of  the  best  expected  value  one  could  achieve  given  advanced  information 
about  the  random  variables.  We  write  the  expected  value  of  this  solution  as 
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J\,  where, 


*1  =  Efrr....,€a[*i(£r,  •••,  6)]»  (3) 

the  expectation,  “E”,  is  with  respect  to  the  random  variables,  f2, . . fr,  and 
Zi  is  defined  as  in  DLP  with  parameters  f2,  •  •  •,  fr- 

An  alternative  and  more  realistic  approach  to  the  stochastic  problem  is 
to  consider  that  during  each  period  t  the  value  ft  is  known  and  a  decision 
on  zt  must  be  made  without  knowledge  of  the  realizations  of  the  future 
periods’  uncertainties.  This  is  known  as  the  “here-and-now”  solution  because 
it  reflects  the  need  for  current  decisions.  The  program  can  be  written  as 

min  z2  =  +  E4a[e2i2  -f  E€,[c3i3  ++  EfT[crzr]] 


subject  to  Aii i  =  &i, 

—Sill  +A2Z2  =  £2, 

—Bt—iZt—i  -f-Ar*r  —  fr. 

it  >0, 

6  £  Bt» 

for  t  =  2, . . . ,  T. 

(4) 


Problem  (4)  states  that  for  t  =  T,  $r— 1  is  given  ,  fr  =  (T  is  observed 
and  >  0  chosen  so  that  ct*t  is  minimum.  Assuming  that  ir  will  be  so 
chosen,  4r— 1  is  chosen  so  that,  cr— iir— 1  +  EfT[crir]  is  minimum  given 
4r— 2  and  fr— 1  observed,  etc. 

In  general,  the  decision  process  for  an  actual  implementation  proceeds 
as  follows: 

A.  A  decision,  2i,  is  made  and  implemented. 

A 

B.  A  realization,  f2,  of  f2  is  observed,  and  a  decision,  t2,  is  made  and 
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T 


1  ■»« 

. .  — 


implemented. 

C.  The  process  is  repeated  to  find  each  t\  in  period  t,  given  the  past  decision 
xt—i  and  the  outcome  of  the  random  vector  fr— i- 

We  are  going  to  investigate  the  effects  of  using  different  methods  for 
determining  the  decisions,  it.  Not  all  of  these  methods  exactly  solve  (4),  so 
we  must  evaluate  the  expected  value  of  the  objective  function  for  solutions 
by  each  method.  For  a  given  method,  n,  of  choosing  xlt  42,  ■  ■  ■,  %t,  we  first 
define  a  function  of  the  random  variables,  ft,  which  gives  the  cost  of  using 
those  decisions.  We  write  this  as 

*($,£  I  **)  =  Cl  •  Xi(fl) 

+  c2  •  £a(£i,  £2  I  /*)  H - h  ct  ■  *r(£i,  •  •  -,  *r— 1,  £2,  •  •  •,  £r— 1  I  m), 

(5) 

where  it(ii,...,  it—  1,  f2,  •  •  •>  £.*  I  V1)  is  the  decision  chosen  by  n  given  the 
previous  decisions,  i lt . . .,  it— 1,  and  observations,  f2, . . .,  ft. 

We  can  then  take  expectations  over  the  random  vectors  to  determine  the 
expected  cost  of  the  solution  found  by  method  We  write  this  as 

*0 <)=  E*M*,ei/0].  (6) 

An  exact  (but  expensive)  method  for  finding  an  optimal  solution  to  (4) 
is  to  proceed  by  a  dynamic  programming  scheme  with  backward  iteration. 
We  will  call  this  “Method  2”.  We  begin  by  setting  the  terminal  valuation 
function: 


zf(xT-  1(fr)=  mincj-ir 

subject  to  Aj'Xj'=  fr  ~f"  Bt—  i*r— 1,  (7) 


0 
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and  let  z\ (zr—  1)  =  E^txJ (xt— i,  £r)]-  We  further  define  the  recursion  on 
the  valuation  function  as: 

4(*t-i,£t)=  min  ctxt  +  z*2+1(xt) 

subject  to  Atxt=  £ t  +  (8) 

xt>  0 

and  let  zfat—  1)  =  E{«[*2(zf— i>  ft)]-  We  finally  arrive  at 
ar£(6i)  =  min  C1X1  -f-  z[x\ ) 

subject  to  Aiii=  bi,  (9) 

h>  0. 

Method  /i  =  2  leads  to  a  sequence  of  decisions,  xt{n  =  2).  We  can  take 
expectations  of  the  outcome  of  these  decisions  as  in  (6).  This  yields 

22  =  E t[z2{x,£  |/i=2)].  (10) 

We  observe  that  for  some  outcome  £  =  ( £2l . . .,  £T), 

2a(*.f  I  M  =  2)  =  cMzl)  +  c2x2^(xi,^2))  H - h  erM*I(4r-i*$r)) 

■where,  for  each  t,  x^z^it—i,  £t))  is  the  optimal  solution  of  (8).  Hence,  by 
integrating,  we  obtain 

z2  =  z\{b  i). 

So,  for  Method  2,  the  result  in  (9)  is  the  minimum  expected  cost  found  by 
the  method. 

The  expected  value,  z2,  is  the  best  possible  solution  for  situations  in 
which  decisions  must  be  implemented  over  time.  (For  details  on  this  result, 


see  Chapter  VI  on  dynamic  programming.)  Enumerating  the  states,  zt— 1, 
in  general,  can  be  very  difficult,  especially  when  £*  can  have  a  continuous 
distribution.  For  this  reason,  we  assume  that  either  the  distribution  of  &  is 
discrete  for  all  t  or  that  we  can  approximate  the  continuous  distribution  by 
a  discrete  sample  of  size  kt  where  kt  is  not  too  large.  We  assume,  therefore, 
that  for  all  t  and  some  E 


p(f.  =  h)  = 


Pt, 

if 

e*  =  e; 

Pt> 

if 

Uy 

II 

p?. 

if 

My 

II 

My 

(ii) 


LO,  otherwise, 

where  by  our  independence  assumptions  the  probabilities  p\  do  not  depend 
on  earlier  outcomes. 

We  have  then  kt  possible  outcomes  for  the  random  right-hand  sides  in 
period  t.  The  outcomes  form  a  tree  of  possible  values  (see  Figure  1.). 


Figure  1.  The  outcome  tree. 
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We  also  define  kt  as  the  total  number  of  possible  outcomes  from  period 
1  to  period  i,  thus 


kt  =  n*=1itr.  (12) 

The  tree  of  outcomes  includes  kt  nodes  in  each  period  t.  We  call  each  node  a 
scenario.  In  period  t  then,  a  scenario  corresponds  to  a  realization  of  outcomes, 

6  =  ii,.  =  h- 

Using  this  framework,  a  descendant  of  a  scenario  j  in  period  t,  is  defined 
as  any  node  in  periods  t  -f- 1  to  T  on  the  branch  connected  to  node  /.  We 
adopt  the  notation  j  for  descendants  of  j.  An  ancestor  of  j  is  then  a  node 
on  the  same  branch  as  j  in  periods  1  to  t  —  1.  We  denote  an  ancestor  of  j 

A 

as  j.  These  definitions  will  apply  in  this  and  all  subsequent  chapters. 

Given  the  discrete  distribution,  we  can  write  (4)  as  an  explicit  linear 
program.  (4)  becomes 


z3  =  min  ciii  +  'EjLiPi^i 
subject  to 

Aix  i 

— #1*1  -\-A2X32 


H - h  E)liPrcrIr 


=  h, 

=  &  for  all  j, 


Bt — 1*3" j  -(- Ay xip  —  for  all  j, 

Xt  >0  for  all  t. 

(SDLP) 

This  program,  SDLP,  is  the  primary  focus  of  our  presentation.  It  repre¬ 
sents  a  formulation  of  the  general  stochastic  dynamic  linear  program  given  a 
discrete  distribution.  By  defining  p\  as  the  probability  of  a  node  in  the  tree  of 
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outcomes,  we  can  also  incorporate  interdependence  of  the  time  periods  into 
the  model. 

SDLP  will  be  analyzed  as  a  structured  linear  program  (presented  in 
Figure  2).  Its  structure  resembles  the  staircase  structure  of  deterministic 
dynamic  linear  programs,  but  the  repetition  of  the  —  B*  blocks  for  each 
descendant  scenario  forms  spikes  below  the  diagonal.  This  property  makes 
the  strict  application  of  staircase  approaches  difficult. 

Another  complication  of  the  stochastic  model  is  that  the  number  of 
blocks,  non-zero  partitions  of  the  coefficient  matrix,  grows  exponentially  with 
the  number  of  periods,  as  we  see  in  (12).  We  will  present  methods  for  solving 
SDLP  that  reduce  the  effects  of  this  complication. 

The  decision  process  of  solving  SDLP  will  be  called  Method  p  =  3.  The 
expected  cost  found  by  this  method  is 

*  =  E[*s(*,£|p  =  3)],  (13) 


Figure  2.  Alternative  structures. 
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An  extreme  simplification  in  SDLP  would  be  to  consider  only  one  pos¬ 
sible  outcome  for  each  random  variable.  A  reasonable  choice  would  be  the 
expectation,  £{.  The  corresponding  program  has  the  same  form  as  DLP  in 
which  the  means,  £ 2 , . . .,  £ T ,  appear  on  the  right-hand  side.  The  solution  of 
this  problem  eliminates  the  increasing  size  problem  of  SDLP  and  is,  therefore, 
the  form  usually  chosen  in  practice.  In  some  instances,  as  we  show  below,  it 
can  even  yield  the  optimal  solution  to  the  stochastic  program. 

We  call  this  deterministic  approach  of  substituting  the  expectations  for 
the  random  variables,  Method  p  =  4.  We  compute  the  expected  cost  of 
Method  4  then  as 

Z4=  Ef[*4(*,£l  A»  =  4)].  (14) 

The  above  four  approaches  to  the  stochastic  dynamic  linear  program  are 
all  related  to  one  another  as  the  following  lemma  states. 

Lemma  1.  The  expected  cost  of  the  approaches  presented  above  are  ordered 

by 

z\  <  *2  <  ?3,  ?1  <  <  *4,  (15) 

and,  for  (discrete)  distributions  without  approximations  z\  <  Zg  =  2a  <  24. 

Proof.  Each  method  can  be  shown  to  improve  upon  the  solution  by  the 
previous  method.  For  <  ?2,  we  observe  that  Z\(x,  £  |  p  =  1)  <  Zi(x,{  j 
p  =  2)  since  Method  1  chose  the  optimum  solution  for  each  £  whereas  p  —  1 
does  not  in  general.  Integration  preserves  the  inequality,  hence,  z\  <  J2- 

The  next  inequality  involves  either  p  =  3,  where  an  incorrect  discrete 
distribution  can  be  used  to  approximate  the  correct  one,  or  p  —  4, where 
the  distribution  is  replaced  by  one  calculated  at  the  expected  value.  Since 
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Method  2  optimises  (8)  for  all  ft,  E ^[zfat-i,  6)]  <  6-i)j 

for  all  t.  Hence,  Z2  <  ^3  or  Z2  <  follows.  | 

The  lemma  includes  the  obvious  result  that,  for  the  expected  value  of 
perfect  information,  EVPI,  where  EVPI=  ?i  —  Z2  if  the  correct  distribution 
is  used  or  ~Z\  —  Z3  if  not, 


EVPI  >  0. 

This  difference  represents  the  maximum  amount  that  one  would  pay  for 
information  about  the  future.  When  the  EVPI  is  low  there  may  be  little 
necessity  in  refining  forecasts,  but,  when  it  is  high,  incomplete  information 
about  the  uncertainties  may  be  costly.  In  the  following  chapter,  we  present 
an  example  of  this  possibility. 

A  second  quantity  that  we  want  to  examine  here  is  the  difference  Z2  — 
Z4  or  Z3  —  Z4,  which  we  call  the  va/ue  of  the  stochastic  solution  ,  VSS. 
VSS  measures  the  benefit  from  solving  a  stochastic  program  over  solving  its 
deterministic  approximation.  A  low  VSS  indicates  that  the  more  complicated 
SDLP  might  not  be  worth  the  extra  effort.  VSS  can  be  bound  without  solving 
SDLP,  as  we  show  in  the  next  chapter. 

We  note  that  in  Lemma  1  to  guarantee  that  Z3  <  Z4,  we  had  to  assume 
that  the  discrete  distribution  was  correct.  If  the  distribution  used  was  only  an 
approximation,  then  it  is  possible  to  make  an  estimate  of  the  distribution  that 
would  lead  to  Z3  >  Z4.  This  anomaly  can  occur  because  some  scenarios  could 
lie  under  sections  of  the  piecewise  linear  curve,  z^xt—  i,  ft),  that  lead  to  high 
penalties.  We  discuss  these  scenario  results  more  closely  in  Chapter  2.  For 
our  purposes,  we  assume  that  the  discrete  distribution  used  suffices  because 
more  information  about  the  distribution  is  not  available.  The  consideration 
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of  individual  scenarios  is  the  only  alternative. 


Given  these  assumptions  about  the  distribution,  we  would  still  like  to 
know  when  SDLP  should  be  solved.  To  show  when  SDLP  need  not  be  solved, 
in  the  next  section,  we  present  some  conditions  that  imply  VSS=0.  In  the 
following  chapter,  we  will  give  bounds  on  VSS  that  also  aid  in  evaluating 
whether  SDLP  is  worth  solving. 

3.  Optimal  Deterministic  Solutions 

When  the  numerical  costs  of  solving  a  stochastic  problem  are  high,  a 
deterministic  solution  technique  is  attractive.  Since  decisions  often  cannot 
wait  for  the  detailed  analysis  of  all  future  possibilities,  the  method  based 
on  assuming  some  "best  guess”  of  the  future  environment,  is  most  often  the 
one  implemented.  In  fact,  the  even  simpler  policy  of  using  a  myopic  solution 
may  provide  a  good  basis  for  decisions.  By  finding  conditions  for  VSS— 0, 
we  are  trying  to  avoid  the  effort  of  correctly  solving  a  general  stochastic 
program.  With  the  conditions  below,  we  can  check  whether  the  stochastic 
program  need  be  tried  at  all.  The  following  lemma,  a  well-known  result  from 
sensitivity  analysis,  is  fundamental  in  our  development. 

Lemma  2.  Let  B  be  an  optimal  basis  for  DLP  with  £  =  (£a,  £s> . . .,  £r).  If 
B  remains  feasible  for  all  £  6  H,  then  B  is  an  optimal  basis  for  DLP  for  all 

t- 

Proof.  Partition  the  coefficient  matrix  and  cost  row  according  to  basic 
variables,  xb,  and  non-basic  variables,  Xff.  DLP  becomes 
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min  ^bzb~\"cszn 
subject  to  Bxb-\-Nxn  = 

ZB  >  0, 

ZN  >  0, 


(16) 


where  we  use  the  notation  t/  to  indicate  the  transpose  of  v,  so  (6j,  £)*  is  the 
column  vector  of  right-hand  sides  in  DLP. 

For  B  optimal,  there  exist  prices,  jr,  such  that 


*B  =  cBl 
icN  <  Cff, 


(H) 


and 


xB  =  B-Hh,DT  >o, 

xn  >  0. 


(18) 


If  xB  remains  feasible  for  all  £  in  (18),  (17)  still  holds,  guaranteeing  dual 
feasibility  and  complementarity.  Hence,  B  is  still  an  optimal  basis.  | 


The  next  problem  we  might  encounter  is  that  of  testing  whether  B  is 
indeed  feasible  for  all  values  of  £.  An  enumeration  of  all  possible  £  is  not 
necessary.  Garstka  and  Rutenberg  (25]  showed  that  simple  computations  for 
many  practical  problems,  could  be  performed  quickly  to  find  the  probability 
that  a  given  basis  is  optimal.  Their  process  involves  fixing  some  components 
in  the  lattice  of  discrete  values  of  £  and  then  finding  the  feasible  range  for  the 
remaining  components.  This  method  also  proves  valuable  in  the  subproblem 
solutions  we  investigate  in  Chapter  3. 

To  use  a  basis  which  satisfies  the  conditions  in  Lemma  2  in  SDLP,  still 
other  conditions  must  be  met.  The  next  lemma  helps  us  find  these  conditions. 
For  this  lemma,  we  will  use  a  solution  from  DLP  in  SDLP.  We  do  this  by 
letting  the  set  of  basic  activities  in  DLP,  {xf  :  t  =  1, . . .,  T},  be  repeated  to 
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form  a  basic  set  in  SDLP.  This  basic  set  in  SDLP  is  {xf1  :  j  =  1,  ...,kt;  t  = 
1  where  xf1  =  xf  for  all  j. 

Lemma  3.  Let  the  set  of  activities  for  a  feasible  basis,  B,  in  DLP  be 
{if,..., if)  where  each  xf  represents  activities  from  period  t.  Also,  let 
SDLP  have  at  least  two  distinct  new  scenarios  at  each  period  (i.e.,  kt  >  2). 
The  activities,  { xf '*},  where  xf 1  —  xf  for  all  j,  form  a  feasible  basis  in 
SDLP  if  and  only  if  xf  consists  of  m(t)  activities  for  all  t. 

Proof.  (See  Figure  3.)  Essentially  we  shall  show  that  if  the  count  on  the 
number  of  basic  xt  is  not  m(t)  for  all  t  that  the  corresponding  SDLP  candidate 
for  basis  will  be  singular.  For  the  necessity  of  the  condition,  first  let  p[t)  be 
the  number  of  elements  in  xf.  Assume  x(t)  does  not  have  m{t)  elements  for 
all  t.  Thus,  there  exists  some  p(t)  >  m[t).  (  If  not,  since  m(t)  =  m, 

p{t)  =  m(t)  for  all  t.) 

We  note  that,  for  p(t)  >  m(t),  t  <  T.  This  is  true  because,  if  fi(T)  > 
m(T),  then  m(t),  which  contradicts  the  fact  that  {xf1} 

corresponds  to  a  basis.  Set  t?  —  min  { t :  p{t)  >  m(£)}.  We  note  that  there 
exists  no  tf1  <  If  such  that  p(t,/)  <  Again,  this  would  mean  the  basis 

was  not  of  full  rank.  Now,  we  have  ]£t=i  /»(*)  ~  Z)!=i  m(t)  =  be  ’  St,  where 
St  =  lift)  ~  Ht1)-  But,  for  t  >  P,  £f=t,+1  M(0  —  m(t)  =  ~h>  • 

^?+i  •  St,  since  each  deficiency  is  repeated  £*'+1  times.  Hence,  A*(0  — 
Xf<(l — ift'+i)  n*(t),  which,  by  our  assumption,  implies  the  columns 

of  the  activities,  {xf},  do  not  form  a  basis  in  SDLP.  To  show  sufficiency,  first 
note  that  if  xf  has  m{t)  elements  in  each  period  t,  then,  for  all  t,  there  exists 
a  square  non-singular  partition  of  the  basis,  Btl  with  columns  and  rows  only 
in  t.  (If  not,  B  does  not  span  the  row  space  in  period  t.)  Therefore,  in  SDLP, 
the  set  of  columns,  {Bt[j}},  is  linearly  independent.  By  construction,  the  xf 
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Wl 

correspond  to  m(l)  -f  z2t=a  ktm(t)  columns,  so  the  activities  form  a  basis. 

I 


From  this  lemma,  we  obtain  our  result  as  stated  in  the  following  theorem. 

Theorem.  If  the  optimal  basis,  B,  for  the  program  DLP,  with  (  =  J,  is 
feasible  for  all  (  £%,  and  if  B  has  as  many  columns  as  there  are  rows  in  each 
period,  then  the  set  of  activities  in  B  forms  an  optimal  basis  in  SDLP,  and 
VSS  =  0. 

Proof.  First,  to  show  primal  feasibility,  let  Af*  be  the  square  non-singular 
submatrix  associated  with  the  activities  xf  in  scenario  j.  For  all  t  and  j,  we 
have 


tf'=(A*')  ‘(if  +  ,),  (19) 

which  is  the  same  value  as  if  in  DLP  for  £*  =  Hence,  ify  >  0  by 
Lemma  1  for  all  j  and  t. 

Let  wf  be  the  dual  variables  in  DLP  for  the  basis,  B.  Next,  define 
irf1  —jPvf.  In  period  T,  we  obtain 


□  □ 

□  □ 

□  a.  □ 


Figure  S.  Example  of  (a)  a  basis  and  (b)  non-spanning  columns  for  a  2  period  example. 
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i 

*=*■ 


«f' = A? 


(20) 


and 

*pA“‘  <  pic?,  (21) 

•where  the  c?  are  the  non-basic  costs  in  scenario  j  for  j  =  1,..., kj>.  In 
general,  for  jr f*  in  SDLP,  we  have 


-  £  *f+iB?y  +  =  ~ P^f+i#?  +  P*xf  A?  =  tfcf  (22) 

and 

-  E  +  'f'A"’  <  cf-  (23) 

7-i 

By  (22)  and  (23),  for  zf*  feasible,  B  is  an  optimal  basis  in  SDLP.  (19) 
implies  that  the  solution  in  SDLP  given  scenario  j  is  the  same  as  in  DLP. 
Hence,  23  =  ?4,  and  VSS  =  0.  We  further  note  that  this  also  implies  the 
optimal  solution  is  myopic.  | 

This  result  gives  us  a  method  to  check  for  deterministic  optimality,  but 
it  may  be  difficult  to  satisfy  these  conditions  in  practical  examples.  Even 
when  not  satisfied,  they  could  be  useful,  however,  in  finding  an  optimal 
solution  with  confidence  a,  where  a  is  the  probability  of  being  optimal.  The 
Garstka  and  Hutenberg  procedure  mentioned  above  would  be  useful  for  these 
computations.  (This  could  also  be  applied  to  a  problem  of  the  form  of  a 
chance-constrained  program  as  in  Charnes  and  Cooper  [14],  but  we  wish  to 
restrict  our  development  to  the  recourse  problem,  SDLP. ) 


I? 
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The  next  result  presents  an  alternative  set  of  conditions  that  may  prove 
useful  when  the  conditions  in  our  theorem  are  not  satisfied.  We  state  them 


as  a  corollary. 

Corollary.  Let  {5(t)}  be  a  family  of  bases  for  DLP,  where  B(i)  is  optimal  for 
£«  €  S.  Assume  also  that,  for  all  periods  t  and  nodes  j  of  scenario  i,  the  set 
of  basic  activities,  ,  is  the  same  for  all  i(j)  that  include  the  nodes, 

{J,  :  l  —  1,  . .  the  descendants  of  j.  If  each  B(i)  additionally  has  square 
blocks  in  each  period,  then  the  set  of  basic  activities  chosen  from  (B(t)}  is 
optimal  in  SDLP,  and  VSS  =  0. 

Proof.  Since  is  the  same  for  all  j,  we  can  define  a  set  of  activities  for 

SDLP  as  xf 1  —  zf(*W  for  all  t  and  j.  Now,  for  primal  feasibility,  we  again 
have  (19)  for  all  t  and  j,  so  xfi  >  0. 

For  the  dual,  define 


Hence,  at  period  T, 


(24) 


irpAf*  =  p*cf, 
and 

ir^/A^  <  p*c]p. 

For  general  t,  we  have 


(25) 

(26) 


-  E  + *?'A?>  =  E  rt— SF* + 


7=1 


(27) 
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and 


Ej+i 

£  * ?l,B?  +  xf'A?'  <  pic?.  (28) 

7=i 

(25),  (26),  (27),  and  (29)  give  us  dual  feasibility  and  complementarity, 
proving  that  the  set  of  variables  {xf*}  is  optimal.  Again,  from  (19),  the 
values  are  the  same  as  in  solving  any  deterministic  form  DLP,  so  VSS  =  0. 

I 

The  corollary  gives  us  more  conditions  for  finding  the  optimal  solution 
to  SDLP  without  actually  solving  it.  An  example  of  a  model  which  meets 
these  requirements  is  the  Hotelling-Nordhaus  model  of  exhaustible  resources 
(see  [34]  and  [46])  and  its  extension  by  Chao  [12]. 

In  Chao’s  model,  a  dynamic  production  schedule  is  chosen  to  minimize 
the  cost  of  satisfying  an  increasing  sequence  of  demand  requirements  over 
time.  The  demands  may  be  satisfied  by  any  of  m —  1  technologies,  each  using 
one  distinct  resource,  with  finite  availability  and  one  “backstop”  technology 
with  no  resource  limits.  The  program  is 

m  oo  m  T 

min  ]P  fiUiyu  +  £  ]C  (29) 

»=i  t=o  «=i  t=o 

subject  to 

Vil  —  R ii  *  1,2,...,  T71, 

i  y»t=  TP*,  t  —  1,2 ,  — ,  T, 

v*.h-i=  V*  +  ~ 

Vit>  0,t  =  0,1,..., 

Zit>  0,t  =  1,2, 
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where  y,t  is  the  amount  of  the  demand,  Dt,  satisfied  by  resource  t  at  time  t, 
Xu  is  the  amount  of  resource  i  committed  at  t  so  it  may  be  extracted  later, 
Ci  is  the  current  cost  of  technology  t,  Jfc,  is  the  capital  cost  of  t,  (3  is  the 
discount  factor,  tf*  is  the  extraction  rate,  and  R,  is  the  initial  availability  of 
the  resource  used  in  technology  t. 

Chao  showed  that,  for  this  model,  a  myopic  solution  is  optimal  for  all 
future  demands  and  supplies.  This  solution  implies  that  a  family  of  bases, 
exists  that  satisfies  the  conditions  of  the  corollary.  Therefore,  the 
stochastic  program  for  (29),  in  which,  Dt  and  Ri  are  random,  has  a  deter¬ 
ministic  and  myopic  solution  and  VSS  =  0.  We  note  that  this  very  simple 
model  can  be  modified  so  that  VSS  grows.  Chao  explored  the  case  of  price- 
responsive  demands  and  found  that,  with  a  “sufficiently  high”  discount  rate, 
the  optimal  decisions  are  still  insensitive  to  “distant-future”  uncertainties. 
Our  example  in  Chapter  II  shows  how  near  future  uncertainties  can  greatly 
affect  current  decisions,  also  making  VSS  high. 

This  chapter  has  described  the  value  of  information  in  the  consideration 
of  decisions  made  over  time.  We  presented  the  program,  SDLP  as  a  method 
for  incorporating  uncertainties  into  a  decision  process  and  explored  the  pos¬ 
sibilities  for  finding  a  solution  to  SDLP  without  solving  the  full  problem. 
The  value  of  the  stochastic  solution  is,  however,  not  always  low. 
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Chapter  II 


The  Nature  of  the  Stochastic  Solution 


1.  Introduction 

In  Chapter  I,  we  presented  conditions,  under  which,  the  optimal  solution 
to  a  certain  deterministic  program  is  the  solution  to  the  stochastic  dynamic 
linear  program.  The  deterministic  model  used  the  expected  values  of  the 
random  right-hand  sides.  .  Unfortunately,  the  great  majority  of  stochastic 
problems  do  not  meet  the  certainty  equivalence  criterion,  ie.  the  sufficient 
conditions  for  deterministic  optimality.  When  a  model  fails  the  conditions 
of  Chapter  I,  it  would  be  desirable  to  solve  the  stochastic  problem  directly. 
In  this  chapter,  we  explore  the  value  of  that  solution  and  the  costs  that  can 
arise  from  not  finding  the  optimal  stochastic  solution. 

We  shall  concentrate  on  decisions  based  upon  “risk  neutrality”,  meaning 
we  wish  to  optimize  the  expected  value  of  our  policy  decisions.  Alternatively, 
the  decision  maker  might  want  to  minimize  the  probability  of  a  catastrophic 
loss  or,  otherwise,  reduce  the  variance  of  his  expected  utility.  These  at¬ 
tributes  could  be  reflected  in  a  carefully  defined  nonlinear  utility  function  or 
in  penalties  placed  on  the  less  attractive  scenarios.  In  this  discussion,  we  do 
not  consider  such  specifications.  Our  consideration  of  linear  models  should 
be,  however,  sufficiently  general  to  allow  for  further  analysis  in  this  area. 

This  chapter  begins  with  a  discussion  of  bounds  on  the  expected  value 
of  perfect  information.  It  then  proceeds  in  Section  3  to  examine  uses  of  the 
deterministic  optimization  of  different  scenarios.  In  Section  4,  the  possibilities 
for  combining  these  solutions  and  the  inherent  difficulties  in  the  stochastic 
program  are  discussed.  Sections  5  then  presents  examples  of  these  problems. 
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The  chapter  concludes  in  Section  6  with  a  suggestion  for  a  general  strategy 
to  be  applied  in  linear  optimisation  under  uncertainty. 

2.  Bounds  on  the  Expeeted  Value  of  Perfeet  Information 

We  discussed  above  the  expected  value  of  perfect  information  (EVP I) 
and  value  of  the  stochastic  solution  (VSS)  and  showed  examples  when  these 
quantities  might  be  zero.  When  the  conditions  for  deterministic  optimality 
are  not  met,  we  would  like  to  have  simple  bounds  on  the  EVPI  that  may 
help  us  determine  the  worth  of  solving  the  stochastic  program.  If  the  EVPI 
and  VSS  are  bounded  within  a  tight  range,  it  may  be  adequate  to  use  a 
deterministic  approach  to  the  problem  instead  of  following  an  expensive 
stochastic  method. 

The  expected  value,  Z,  of  the  objective  function,  2(f),  can  be  bounded 
because  of  its  convexity.  Madansky  (42]  and  later  Huang,  Ziemba  and  Ben- 
Tal  [35]  examined  this  property  using  the  theory  of  moment  spaces  to  bound 
the  expectation.  Their  work  rests  on  the  following  result. 

Lemma  1.  The  objective  function,  z($),  in  (1.3)  is  a  convex  and  continuous 
function  of  £,  the  right-hand  side. 

Proof.  See  Madansky  [42].  | 

This  result  then  allows  the  application  of  Jensen’s  inequality  for  convex 
functions.  Directly  from  this,  we  have 

=  E [z(f)]  >  z(  E(0)  =  *4,  (1) 

giving  a  lower  bound  on  the  perfect  information  solution,  Zj.  (In  this  analysis, 
we  use  the  definitions  of  z\,  zg,  zg,  and  Z4  from  Chapter  I.) 

An  upper  bound  also  can  be  found  for  Zj.  We  present  here  only  the  one 
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dimensional  case.  The  multidimensional  case  can  be  found  in  Madansky  [42]. 
The  bound  is  known  as  the  Edm undson-Madansky  inequality,  and  it  states 
that,  for  £  €  [a,  6]  C  Ut\/x  =  E[£], 


=  EWf)]  <  )*W  +  W 

i  i  i  p  *j  •  i  rtfpl  mi  *  •  i  m  •  ft  i  ^ 

wnere  oy  aenmuon  nere  fx  =  ims  is  snown  easily,  since  ior  any  1 1 

[a,  6],  t  =  ( a(b  —  t)  -f-  (t  —  a)b)f{b  —  a)  =  Xo  +  \b,  and,  for  z{t)  convex, 
z[t)  <  \z(a)  -f-  \z(b),  which  by  integrating  yields  (2).  We  can  now  state  (1) 
and  (2)  as 

Theorem  1.  For  £  £  [a,  6]  C  JR1,  fx  =  E(£),  and  z  as  defined  above, 

(6  —  n)/{b  —  a)z(a)  +  (/i  —  o)/(5  —  d)z{b)  >  Zi  >  z[fx). 

Huang,  Ziemba  and  Ben-Tal  carry  these  principles  further.  They  sub¬ 
divide  the  interval  [o,  6]  and  apply  successively  finer  approximations,  which 
they  show  approach  the  expectation.  This  method  makes  possible  the  refinement 
of  the  EVPI  to  whatever  level  is  desired. 

Another  method  for  computing  bounds  on  the  EVPI  was  presented  by 
Avriel  and  Williams  [5].  They  showed  that 


0  <  EVPI  <  z2  —  z4  <  Z4  —  z4,  (3) 

where  Z2  is  the  best  stochastic  solution  and  z4  is  the  expected  value  deter¬ 
ministic  solution,  as  in  Chapter  I.  They  also  show  that,  for  z(()  differentiable, 
the  bound  using  the  expected  value  of  £  (as  in  z4)  is  the  tightest  possible. 

These  approaches  give  us  methods  for  estimates  of  the  benefit  we  may 
gain  from  knowledge  about  the  random  variables.  We  concern  ourselves  here 


with  the  VSS,  the  value  of  solving  a  stochastic  program  over  solving  a  deter¬ 
ministic  one.  The  EVPI  is  used  to  find  the  value  of  additional  information, 
but,  in  looking  at  the  VSS,  we  assume  that  no  more  information  is  available. 
We  ask:  what  given  our  present  state  of  knowledge,  is  the  value  of  solving  a 
large  stochastic  program? 

Similar  results  to  those  above  can  be  found  for  VSS.  We  can  bound  VSS 
as  in  the  following  theorem. 

Theorem  J.  The  value  of  the  stochastic  solution  (VSS  =  Z4  —  ?3,  as  defined 
in  Chapter  I),  satisfies  the  inequalities 

0  <  VSS  <  Z4  ~  *4-  (4) 

Proof.  We  showed  VSS  >  0  in  Chapter  I.  It  suffices  to  show  z4  <  2S-  We 
had  Zi  <  <  ?3-  Now,  z4  —  z(  E(£))  and  =  E(z(£)),  so,  for  z(£) 

convex,  by  Jensen’s  inequality,  z4  <  .  The  result  follows.  | 

This  bound  can  prove  useful  in  estimating  the  benefit  of  the  stochastic 
program,  but,  if  it  remains  high,  further  analysis  may  be  required.  As  a 
first  step  in  solving  the  stochastic  program,  we  may  find  other  deterministic 
solutions  corresponding  to  different  scenarios  or  outcomes  of  the  random 
variables.  We  describe  this  approach  in  the  next  section. 

3.  The  Scenario  Approach 

In  evaluating  the  perfect  information  or  “wait-and-see"  solution,  ~Z\,  a 
solution  to  the  program,  zj({),  must  be  found  for  each  possible  outcome  of 
the  random  vector,  £.  The  scenario  approach,  also  known  as  “modified  wait- 
and-see"  in  Gunderson,  Morris,  and  Thompson  [31],  involves  solving  several 
of  these  deterministic  programs,  evaluating  the  expected  cost  of  using  the 


strategy  that  is  optimal  for  each  scenario,  and  choosing  the  strategy  that 
minimizes  this  expected  cost.  In  many  cases,  one  of  the  first  period  bases 
dominates  the  others,  and  the  choice  for  a  decision  is  clear.  As  Gunderson,  et 
al,  emphasize,  this  method  can  be  quite  responsive  to  management  concerns 
and  may  prove  very  useful  since  it  presents  alternative  possibilities  and  risks 
in  a  compact  and  easily  understandable  form. 

In  the  scenario  approach,  we  first  choose  a  set  of  possible  outcomes  for 
£,  which  we  call,  H,  where 

a  =  (5) 

For  each  £°  6  H,  we  find  the  optimal  solution,  x* ,  for  z(£°).  Then,  we 
compute 

f({°)  =  E[f(**({°),  f)]  (6) 

for  each  £°,  where  f(z*(£°),  £)  is  the  resulting  objective  value  from  using  x* 
when  (  actually  occurs. 

Next,  we  find 

f*  =  min  €oeg  f(£°).  (7) 

This  value  represents  the  least  expected  cost  from  using  the  solution  of  a 
deterministic  program. 

This  implies  that  a  deterministic  problem  other  than  the  expected  value 
problem  may  result  in  a  better  solution.  This  is  because  less  penalty  may 
be  incurred  by  following  a  piecewise  linear  section  of  the  objective  function 
other  than  that  that  covers  the  expected  value.  (See  Figure  1.) 

We  can  also  bound  f*  as  in  the  following  theorem. 


Theorem  S.  For  B  discrete,  the  best  scenario  solution  ,  f*,  satisfies 


*3  <  f*  <  *4 


(8) 


Proof.  We  know 

f-  =  min  E!f(i(f0)),  f]  <  E[f(*(D,  f))  =  J, 
and,  for  a*  optimal  in  C', 

Jg3  =  min  *eB  E[f(*(f)),  e)  <  EM**),  £]  =  f*. 

Hence,  the  result.| 

The  inequalities  in  (8)  show  that  the  scenario  approach  may  be  useful 
in  finding  a  closer  approximation  to  the  solution  of  the  stochastic  program. 
?*  may  be  especially  valuable  when  f*  —  z*  is  small,  since  it  also  bounds 
Zs  —  Z+  and  may  show  that  solving  the  stochastic  program  is  unnecessary. 


z(x(5) ,5) 
z(x(£) ,5) 

z(x(5) ,5) 


Figure  1.  The  expected  value  of  ()  is  greater  than  that  of  z(x(£),  £). 


This  outcome  is  most  likely  when  the  optimal  decisions  in  the  first  period 
(those  that  are  actually  implemented)  correspond  to  the  same  basis  for  most 
of  the  scenarios. 

The  scenario  approach  can  also  be  used  to  find  activities  that  may  be 
basic  in  the  optimal  solution  of  the  stochastic  program.  If  the  set  of  optimal 
basic  activities  remains  fairly  constant  for  the  possible  realizations  of  £,  then 
this  set  of  activities  may  be  optimal  in  the  stochastic  program.  We  discuss 
this  possibility  in  the  next  section. 

4.  Combining  Scenarios 

A  natural  approach  to  solving  the  stochastic  program  would  be  to  use 
the  optimal  solutions  of  the  different  scenarios  and  to  combine  them  in  an 
appropriate  manner.  The  proper  combination  may,  however,  be  quite  difficult 
to  find  and  may  lead  to  as  much  effort  as  solving  the  stochastic  program 
directly.  We  present  below  the  problems  inherent  in  combining  scenarios  and 
some  situations,  in  which,  the  optimal  combination  may  be  found  directly. 

In  this  analysis,  we  restrict  ourselves  to  a  two  period  case,  for  which  there 
are  only  two  scenarios  considered.  The  difficulties  involved  in  this  example 
are  typical  of  all  stochastic  programs,  so  we  present  this  case  because  of 
its  simplicity.  The  results  may  be  easily  generalized  to  more  periods  and 
scenarios. 

We  begin  by  defining  two  scenario  problems  as 


min  ^(f1)  =  CiXi-j-c3X2 
subject  to  A1X1  =  61 

—  BiIi+A2X2=  £2 
*1>*2  >  0, 


(51) 


and 
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(52) 


min  z((2)  =  CiZi-\-C2Z2 
subject  to  AiZi  =  bx 

—BiZi-\-A2Z2—  £2 

Z\t  z2  >  0. 

Solving  SI  and  S2  yields  the  optimal  solutions,  z1’*  and  z2'*,  and  optimal 
dual  prices  (jt1**; ct1**)  and  {x2‘*]02,*)>  'where  z*  =  it*  ■  bx  +  o*  •  £.  We  also 
define  the  following  index  sets: 

yS1  =  {j  :  Ai(*,j)  is  basic  in  SI  } 


and 

/ 3 2  =  { j  :  Ai(*,j)  is  basic  in  S2  }.  (9) 

1  2 

The  complements  of  yfl1  and  ff2  are  defined  as  ]5  and  ,  respectively. 

Now,  our  actual  goal  is  to  solve  the  following  two-period  version  of  SDLP 


min  cxzx  +p1c2x^  +f>2c2*l 
subject  to  AiZi  = 

-Bizi  +Mz\  (SDLP2) 

—B\Z\  +A2z|  =  £§ 

We  would  like  to  use  the  activities  in  fil  and  fi2  as  the  optimal  basis 
of  SDLP2.  In  other  words,  for  /?  =  {j  :  Ai(*,j)  is  basic  in  SDLP2  },  we 
are  looking  for  /?  C  01  |J  fi2.  Unfortunately,  this  is  not  always  possible. 
The  difficulty  results  from  the  properties  of  the  basis  in  SDLP2.  In  order  to 
maintain  full  rank  of  the  basis  in  this  program,  the  optimal  basic  activities 


from  (SI)  and  (S2)  cannot  all  be  in  the  optimal  basis  of  SDLP2,  unless  the 
square  block  situation  we  discussed  in  Chapter  I  holds. 

In  general,  when  we  attempt  to  combine  scenarios,  we  face  a  duality  gap 
as  the  following  lemma  states. 

Lemma  2.  For  z* ,  the  optimal  value  of  SDLP2, 

z  =  {p1o1'*  +pV'*)6+p1jr1'*^-t-p2jr2'‘^l 
<  (10) 
<  Ciip'zl'*  +  p2zl'*)  +  pciiz’*  +  pae2z2'*  =Z, 

where  x\’*  and  z2’*  are  the  optimal  primal  solutions  for  SI  and  S2,  i2’*  and 
ij'*  are  the  values  of  the  second  period  basic  variables  in  SI  and  S2  chosen  to 
satisfy  A2x\'*  =  +  ^(P1*!'*  +p2z?'*)  and  A2z2'*  =  £2  +  SitP1^'*  + 

p2zl'*),  [if  z\‘*  or  z2'*  <  0,  sefZ  =  -fooj,  and  x1'*,  x2'* ,  o1’*,  and  <r2'*  are 
the  optimal  dual  values  for  SI  and  S2.  Furthermore,  the  duality  gap,  if  z\'* 
and  ij’*  are  feasible,  is 


*-*  =  (*!-  +  x2‘*B\)p2plz1'*  +  (c\  —  ox‘*  A2  +  xl,*Bi)pxp2z2,*t 

(11) 

where  cj  =  { cy  :  j  G  P1C\P2}>  c2  =  {ey  :  j  €  pl  DPa},  A\  =  {ay  :  oy 
is  a  column  of  A\  and  j  £  P1  f]P  },  Aj  =  {oy  :  oy  is  a  column  of  A\  and 
j  €  ft1  Cl  P2},  and  and  B\  are  columns  of  B\  corresponding  to  A\  and 
A2,  respectively. 

Proof.  First,  we  show  that  (p1o1»*  p2oa'*;p1x1’*>p2x2'*)  is  a  feasible 

solution  to  the  dual  of  SDLP2.  We  have  jra,*A2  <  e2  and  jt1,*A2  <  c2, 
since  jt1’*  and  ir2,‘  solve  SI  and  S2.  This  also  yields  — xx,*Bi  -|-cr1,*Ai  <  c\ 


-(P1*1'*  +/>V'*)£i  -f  (pV**  +  pV'*)Ai  <  Cj}  (12) 

since  p1,  p2  >  0  and  p1  -j-  p2  =  1.  Therefore,  the  solution  is  dual  feasible, 
hence,  by  duality,  we  hare  z*  >  z. 

Now,  we  consider  (p1  x^+pV'*;  zj**;  x2-*)  and  observe  that  A1(p1i1**+ 
p2i2'*)  =  6,  that,  by  definition,  and  x2'*  satisfy  the  second  set  of  in¬ 
equalities  in  SDLP2,  and,  that,  if  i2**  and  i2'*  are  not  feasible,  2  =  oo. 
Therefore,  if  x2*  >  0  and  i2’*  >  0,  the  solution  is  feasible  and  z*  <  z. 

Therefore,  z  <  z*  <  2.  For  the  expression  of  the  duality  gap,  we 
observe  that 


z  =  (pV-*  +  pV'*)^^*1'*  -I-  p2AiZ2'*) 

+  p1T1-*(A2iJ’*  -p1Biz{'*'  -p2BlZ2'*) 

+  P2jt2-*(A25J*  -p2B!Z2'*) 

=  ^(P1*}’*  +P2*?'*)  +pM(p'zl'*  +  p2x2'*) 

+  p2c2(p1xj**  +  p2*?'*)  -fP^'V^*!'*  +  p2A}z2’*) 

+  p2x2‘*(-p1BjxJ.*  p2BjX2’*)  -j-  p1<71,*(p1A2xJ'*  -f- p2A2x2’*) 

+  A^-p^x}’*  — paB2x2'*)  +  A*zJ-* 

+  p2c|x2'* 

(13) 

where  c°  =  {Cl(j) :  j  e  Z?1  f|0a}.  We  then  have 

z  =  7  -  p2cj(p1xj'*  +  p2x2'*)  -  p^AJ-*  +  p2x2'*) 

-f  (pV-*Al  -  p2x2,*B}  +  pV-‘A2  -  p1?rI-*B2)(p1xi'*  +  A?'*) 
which,  since  e\  •  x}'*  =  0  and  cj  •  x?'*  =  0,  yields 


z  =  2  -  (c[  -  aa-*Aj  +  x2-*B})(p2p1xJ'*) 

-(cj-^Ai  +  ^^yxj.*) 

w*e  r«Mk  is  (11)  folkw*.  | 
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Measuring  the  duality  gap  in  (11)  is  another  way  of  finding  the  value  of 
solving  the  stochastic  problem  directly.  If  the  gap  is  small,  then  we  can  use 
the  simple  weighted  average  of  the  optimal  scenario  values  without  incurring 
a  great  penalty.  A  large  gap  signifies  greater  variance  among  the  scenario 
solutions  and  should  lead  to  a  stochastic  approach.  A  large  difference  can 
also  result  from  infeasibility  of  the  primal  solutions  since  by  definition  in  this 
case,  Z  =  oo.  In  some  instances,  the  gap  may  be  closed  easily  as  we  discuss 
below. 

We  first  consider  the  case  of  /?  =  /31  =  (same  first  period  basic 
activities  in  SI  and  S2),  but  where  |/9|  >  mi,  the  rank  of  Ai.  This  situation, 
which  we  described  in  Chapter  I,  signifies  that  the  same  basic  activities  cannot 
exclusively  be  used  to  solve  SDLP2.  They  will  not  span  the  row  space  of  the 
matrix.  If  we  have  an  inequality  form,  however,  the  weighted  average  of 
optimal  first  period  values  may  be  optimal  in  SDLP2. 

We  write  the  coefficient  matrix  of  SDLP2  as 


A  = 


/  A?  A\  A\  A? 

-B°  -B\  -B\  —B?  A^A%N 
k—b\  —B\  -b\  -b?  a|aJ" 


(16) 


where  A°  =  (Ai(*,;)  :  j  G  PlC\p2),  A\,  A?  are  defined  as  in  Lemma  2, 
A?  =  (A,(.,  j)  :  j  €  PI  f),  Ai,  A!  are  the  basic  second  period  columns 
in  SI  and  S2,  respectively,  the  columns  of  A^N  and  are  non-basic,  and 
the  B)' s  are  defined  correspondingly. 

Now,  we  assume  that  j3\  =  £2  and  that  we  can  replace  the  second 
equality  in  SDLP2  with  an  inequality.  SDLP2  becomes 
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mm 


e?»i  +'!»!  +c"vi  +p1‘  Ivl  +p*4vl  +p'fiNvl,N  +pafiNvlK 


».t.  A°|,°  +A}»J  +Af'»{'  =6, 

-B?»!  -Bill!  -B?y?  +A>» I  +Ai.w„«.w  >  f  1 

-B»»«  -Bis,!  -Bf'„['  +A|y|  +fi"fiN  >  fl 

Vl,  ?2  >  0, 


(SDLP2') 

where  the  variables  (yijife)  replace  the  variables  (21,22)  in  SDLP,  so  that 
we  can  compare  their  values.  First,  we  define  a  solution  to  SDLP2'  by  the 
following.  Let  a  solution  {y\,y\,y\)  be 


This  solution  is  the  same  as  ( x\ z\’ *)•  Next,  define 

i/  =  {*:  2{'*(t)is  basic  in  row  j  of  SI  for  j  >  mi},  (18) 

where  mx  is  the  number  of  rows  in  A\ .  Now,  partition  as 

Ai'v  =  *) :  j  6  *0 

and 

Ai'”  =  (A>W, .) :  J  5?  ■'I-  (19) 

is  non-singular  because  the  basis  in  (17)  is  non-singular. 

We  complete  the  definition  of  y*  by 


fir  =  (A*’1’)— *(#r  +  +  Bi^l'-) 


and 
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y\y  =  ft"  +  +  B\‘vy\'\  (20) 

Here,  y2'*''*  represents  the  slack  variables.  Now,  we  wish  to  show  that  y*  is 
an  optimal  solution  to  SDLP2*.  We  will  have  to  restrict  the  solutions  of  SI 
and  S2  for  this  to  be  true.  We  first  look  at  the  dual  of  SDLP2*. 

Let  jt2'*'  =  (ir2  :  j  E  v)  ==  0  and  x2-F  =  (x2  :  j  gv)  =  P2c2(j4|,1')“1,  so 
we  have  x2i4^  =  Cp2c§;  0)  where  the  objective  row  coefficients  of  the  slacks 
y\v*  are  0. 

Next,  let 


r  l'v = 


—  X 


(21) 


and 


»*•'  =  (tl  -  «K  -  Bfp1ei(^)-,)(Br(A‘'I’)-1A>'-  -  Btr1, 

where  B\  =  (S1(*,  t)  :  i  is  basic  for  row  j  such  that  j  E  v)  and  Bj,  A*  are 
defined  accordingly.  Now,  let  B  =  {B^v A^,v  —  By)~l  and  compute 
er  as 


*=i«r+ w  - 

+  Bf(A‘'I7)-1c‘!(A'’r  +  (Sr  -  BJ(A>'‘')-'^‘')(Ar1A0-‘. 

We  substitute  for  t1>*/  in  (21)  and  find  that,  *old  and  gold,  the  optimal 
dual  prices  in  Si  (or  S2,  since  the  bases  and  objective  functions  are  the  same), 
satisfy 
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(i)  jr2-"  =  0, 

(ii)  *»•>’’  =  p^SlD  + 

(iii)  xv  —  KqLq, 

(iy)  r‘.”  =  p‘(rgLD  -  (p‘/p‘)*tLDA}^4',’r\ 

and  (v)  a  —  tr0L0. 

Hence,  we  have 


-x*rBl  -  ^Bl  -  IT *'*Bt  -  k1,i'Bi  +  crAt 


*oldB i  *oldBi  +  oqldA i  <  ci, 


*oldB  f  *old&  f  +  gold  A*  ~  c*[, 


^old^i  ~  *oldB  i  +  gold  A*  = 


We  also  observe  that 


+  y'.-4-  =  pVS^A^  +  rg^y, 

~"Pa’rOL£Jj42,'/  4“P2’rOL£)^'*' 


=  AL 


and  that,  similarly, 


We  need  only  that 


+  %l'vA^’v  < 
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(28) 


and 

x2,va£,v  -f  fc2,vA^,v  <  p2c\N  (29) 

in  order  for  feasible  y*  to  be  optimal,  since  we  have  already  shown  com¬ 
plementarity  in  (24),  (25),  (26),  and  (27).  (23)  also  shows  dual  feasibility  for 
the  first  period  constraints. 

To  ensure  that  (28)  and  (29)  hold,  we  need  to  have  the  following  equality 
satisfied: 


v  a2,u ( a2,T7\  —  IaN.v  _  u  AI\ 

*OLDa2  i-^2  i  -^2  —  ^OLD-^i 


N,  i/ 


By  the  definition  of  jr1  and  j r2,  (30)  implies  (28),  since 


(30) 


A,y  an,v 


=  p'^OPD^  +  *OLD*?’*) 

+  A*"oU>*S*  - 

<  P '<■?, 


and  (29)  holds  by 


AZ-V  + 


1,v  AN,V 


=  p\*oldA%’V  +  (x'oLD^^r'A?'*)) 
~  P2{xOLD A?'17  +  xOLD^’¥) 


<  P  c{ 


Hence,  we  have  shown  the  following  theorem. 


Theorem  4.  For  y*  as  defined  in  (17)  and  (20),  if  y*  is  a  feasible  solution 
to  SDLP21  and  the  optimal  dual  prices  in  SI  and  S2  are  such  that  (30)  is 
satisfied,  y*  is  an  optimal  solution  to  SDLP2f . 

The  restriction  of  the  prices  in  (30)  guarantees  optimality,  but  (28)  and 
(29)  may  be  true  even  if  (30)  does  not  hold.  In  examining  a  problem  of  this 
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type,  if  (30)  fails,  one  may  want  to  carry  out  the  additional  computations  in 
(28)  and  (29)  before  solving  the  stochastic  problem.  , 

The  difficulty  in  finding  optimality  conditions  for  combining  scenarios 
for  even  a  simple  problem  such  as  SDLP2',  shows  the  importance  of  the 
stochastic  solution.  The  conditions  in  (30)  can  be  generalized  to  allow  for 
penalties  in  satisfying  inequalities  in  the  second  period,  but  the  results  are 
more  restrictive  and  direct  computation  of  the  dual  feasibility  conditions  be¬ 
comes  more  efficient  than  checking  additional  inequalities.  Our  development 
leads  to  the  following  method  for  combining  scenarios. 

(A.)  Combine  the  first  period  scenario  solutions,  x\, . . . ,  zf ,  by  a  simple 
weighted  average,  y\  =  53*=iP<a:i- 

(B.)  Follow  the  branch  of  worst  cases  (what  we  shall  call  the  catastrophe 
branch),  that  is,  the  set  {£*}  where  £ J  =  (£J(1),  £J(2), . . . ,  £J(mt))  such  that 
£}(*)  =  sup  for  t  =  1, 2, . . . ,  TOf.  Using  these  right-hand  side  values, 
we  have  that,  if 


>  £t*, 


then 


for  all  j.  This  procedure  guarantees  primal  feasibility  lor  yj.  We  use  the 
non-singular  blocks  from  this  scenario  to  determine  new  values  based  on  y[ . 

(C.)  Use  similar  blocks  for  the  other  branches  and  maintain  primal 
feasibility  (by,  perhaps,  paying  penalties). 

(D.)  Compute  the  dual  prices  and  check  for  closure  of  the  duality  gap. 
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This  approach  would  be  most  successful  when  the  different  scenarios 
have  nearly  identical  bases.  In  these  cases,  the  noncomplementarity  would 
exist  in  only  a  few  terms.  The  difficulty  of  using  these  optimal  linear  program 
solutions,  however,  is  that  they  all  correspond  to  extreme  point  solutions  and 
may  include  very  different  sets  of  basic  activities.  This  property  could  make 
their  combination  in  a  stochastic  program  most  difficult.  In  the  next  section, 
we  present  small  examples  of  this  occurrence. 

5.  Examples 

The  extreme  point  properties  of  the  basis  in  a  linear  program  are  cru¬ 
cial  in  understanding  stochastic  program  solutions.  Critical  values  of  the 
parameters  limit  the  use  of  different  scenarios.  In  some  cases,  the  implemen¬ 
tation  of  any  deterministic  solution  may  lead  to  heavy  penalties  relative  to 
the  solution  of  the  stochastic  program.  One  example  of  this  occurrence  is 
the  following  linear  program  : 


min  z—  zx 

+4za 

+  E^[  min  yx 

-|-10pa 

|  zx  and  z2] 

subject  to  zx 

-j-*2 

=  1 

—Zi 

+2  za 

+Vi 

=  e 

0<yi 

<2, 

*l,Z2,V2 

>0, 

£  is  Uniform  [0,4]. 

{EX  1) 

We  solve  EX1  for  £  =  (1,3)  and  find  the  expected  value  of  using  the 
optimal  decisions  for  these  scenarios  as 

*T=  E[*j(f)].  (31) 


where 
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(32) 


*?<{)  =  ei*(?)  +  [  min  ,  eay|r(|),  {). 

We  define 


7/  =  7(1)  (33) 

and 

-in  =  *(3).  (34) 

We  want  to  consider  also  the  perfect  information  solution 

Zp  =  E(  [  min  c%  ■  x  +  c2  •  y|£],  (35) 

and  the  stochastic  solution,  7S,  where  we  allow  £  =  1  or  3  with  equal 
probability. 

The  function  z^(£)  for  Z/,  Zu,  z„  and  zp  appears  in  Figure  2.  We  observe 
that  the  optimal  basis  changes  as  £  ranges  over  [0,4]. 

For  £  <  1,  Z\  only  is  in  the  optimal  basic  set  of  variables,  for  1  <  £  <  3, 
{%i,  z2}  is  optimal,  for  2  <  £  <  3,  an  alternative  optimal  set  includes  z2 
alone,  and,  for  3  <  £  <  4,  the  only  optimal  first  period  activity  is  x2.  In 
the  stochastic  solution,  z\  and  z2  must  be  in  the  optimal  basis,  reducing  the 
expected  loss  relative  to  the  perfect  information  solution.  We  find  this  from 
the  figure  as 


7/  —  Zp  =  10.25, 


(36) 


in  —  Zp  —  5.50,  and 
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(37) 


7,  —  7p  =  3.75. 


(38) 


The  losses  in  (36),  (37),  and  (38)  associated  with  this  problem  demonstrate 
the  usefulness  of  the  stochastic  solution.  By  assuming  any  deterministic  value 
for  the  right-hand  side,  a  large  loss  may  result.  The  simple  stochastic  for¬ 
mulation  with  two  possible  £  values,  however,  reduces  the  risk  of  this  situation 
and  enables  us  to  approach  the  perfect  information  solution.  The  stochastic 
solution,  therefore,  lends  resilience  to  the  result.  It  provides  a  rationale  for 
hedging  strategies. 

To  demonstrate  further  the  sensitivity  of  models  to  uncertainty,  we 
return  to  the  exhaustible  resource  model  of  Chapter  I,  (1.29).  We  stated  that 
this  model  had  an  optimal  deterministic  solution,  but,  by  adding  an  uncertain 
return  from  investment  in  exploration,  we  again  arrive  at  a  situation,  in 
which,  every  deterministic  solution  will  be  associated  with  losses  relative  to 


/ 


Figure  2.  Example  costs. 
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the  stochastic  solution. 

The  modified  problem  is 


min  X)<==ic*u»,o  fc»y«,o  + 


subject  to 


+  u{t- 1 


£°  +  “*,o 

< 

Ri,0> 

for  *  =  1,..  .,m, 

]Ct=l  u«',o 

> 

■Do; 

<i  +  <t 

< 

Ri.r> 

for  all  i,  j,  and  t, 

> 

Dt;  for  all  j  and  t, 

— W-1 

= 

r 3  • 
*i ,t> 

for  all  i,j,  and  t, 

t  >  0, >  0,  for  all  i,  j,  and  t, 

( ERM ) 

where  Ritj  Dt,  x{t,  and  u3it  are  as  defined  in  (1.29).  y3it  represents  the 
amount  invested  at  time  t  under  scenario  j,  and  a3it  is  the  return  for  invest¬ 
ment. 

ERM  includes  deterministic  investment,  but  we  can  formulate  an  as¬ 
sociated  stochastic  model  by  allowing  a3'  to  take  on  several  values,  a^’1 , . . .,  aV 
and  restricting  investment  in  yi  t  to  be  only  one  of  the  a3,f  for  each  scenario. 
By  doing  this,  we  maintain  the  random  variable  in  the  right-hand  side  as  a 
constraint  on  y?  t  for  the  different  scenarios. 

By  solving  ERM  for  fixed  values,  the  investment  decision  may  quickly 
swing  from  one  resource  to  another  as  a  different  extreme  point  in  the  linear 
program  becomes  optimal.  The  stochastic  solution  has  more  basic  activities, 
allowing  for  hedging  against  the  different  possible  environments  that  one  may 
face.  To  show  this  property,  we  consider  a  two  period  case  of  ERM  with  three 
alternative  resources.  We  will  call  them  oil,  solar  power,  and  some  other  high 


cost  backstop  technology.  We  use  the  following  inputs  in  Table  1. 


Resources 


Cost(ei)  Availability(Rt) 

Oil 

5 

25 

Solar 

10 

10 

Backstop 

16 

00 

Investment 

Cost  (Jfc,-) 

Return(a,) 

Oil  1 

1 

1 

Oil  2 

1 

0.1 

Solar 

1 

1 

Demand 

Period  1 

15 

Period  2 

25 

Probabilities 

Scenario  1 

0.5 

Scenario  2 

0.5 

Discount  Rate 

P  =  .80 

Table  1.  Two  period  model  inputs. 

The  only  uncertainty  in  this  model  is  on  the  return  for  oil  exploration. 
Investment  in  solar  power  can  be  interpreted  as  the  relatively  certain  amount 
of  capacity  increase  from  investment.  Each  unit  invested  in  oil  in  this  model, 
however,  results  in  either  a  full  unit  increase  or  a  tenth  of  a  unit  increase 
in  oil  availability  according  to  a.  These  two  scenarios  are  assumed  to  occur 
with  equal  probability. 

The  deterministic  models  for  “bad  luck"  (a  =  .1),  “good  luck*  (a  =  1), 
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and  "myopic”  (solving  only  given  the  first  period  availabilities  and  demands) 
were  solved  and  compared  to  the  stochastic  program  solution,  in  which  both 
values  of  a  were  considered  simultaneously.  The  results  were  : 


Model _ Expected  Cost 


Myopic 

275 

Good  Luck 

239 

Bad  Luck 

245 

Stochastic 

231 

Here,  the  stochastic  solution  represents  a  savings  over  the  deterministic 
models  because  its  solution  involved  investment  in  both  oil  and  solar  tech¬ 
nologies,  while  the  deterministic  scenarios  allowed  for  investment  in  only  one. 
The  relative  savings  would  also  increase  with  a  higher  cost  backstop.  It  is 
interesting  to  note  also  that,  with  the  addition  of  investment  in  exploration, 
the  myopic  solution  is  now  far  from  optimal. 

These  examples  have  shown  that  models  can  be  very  sensitive  to  future 
uncertainties.  The  exhaustible  resource  model  demonstrates  the  possibility  of 
sharp  changes  in  decision-making  from  near-term  uncertainties.  We,  there¬ 
fore,  want  to  examine  models  with  a  strategy  that  considers  their  sensitivity 
to  uncertain  parameters.  We  discuss  a  general  method  for  dealing  with  these 
problems  in  the  following  section. 

6.  A  General  Strategy 

Problems  modeled  as  dynamic  linear  programs  often  involve  many  un¬ 
certain  assumptions  about  the  parameters  involved.  In  this  case,  a  solution 
to  the  stochastic  program  is  desired,  but  it  may  be  quite  costly.  We  have 
presented  methods  for  checking  whether  deterministic  solutions  may  be  used. 
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The  strategy  resulting  from  our  development  in  Chapters  1  and  2  follows. 

I.  Solve  the  expected  value  problem  and  check  for  basis  feasibility  and 
a  possible  optimal  deterministic  solution  as  in  Chapter  I. 

II.  If  I  fails,  use  the  properties  presented  in  2.3  to  bound  the  VSS  and 
EVPI. 

III.  If  the  VSS  bound  is  large,  solve  different  deterministic  scenarios 
and  look  for  a  dominant  basis.  Then,  attempt  to  close  the  duality  gap  by 
combining  scenarios. 

IV.  If  the  gap  persists  and  if  no  single  basis  is  indicated,  proceed  to  solve 
SDLP. 

This  procedure  outlines  how  we  would  evaluate  the  worth  of  solving 
successively  more  complicated  problems.  As  our  examples  have  shown,  we 
may  find  that  the  stochastic  program  is  worth  solving,  and  that  its  solution 
may  result  in  a  substantial  savings.  We  might  attempt  to  solve  such  problems 
directly  by  brute-force  methods.  Alternatively,  the  next  three  chapters 
present  methods  that  take  advantage  of  the  structure  of  SDLP  to  reduce 
the  sise  of  these  possibly  very  large  linear  programs. 
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CHAPTER  m 


A  Nested  Decomposition  Algorithm  for  SDLP 


1.  Introduction 

Stochastic  dynamic  linear  programs  have  natural  subdivisions  correspond¬ 
ing  to  the  separate  decisions  made  in  each  period  under  various  scenarios. 
These  divisions  give  SDLP  a  special  structure  that  can  lead  to  improved 
efficiency  in  its  solution.  In  the  next  three  chapters,  -we  present  three  methods 
employing  distinct  optimization  techniques,  each  of  which  exploits  the  pro¬ 
gram’s  structure.  The  structure  is  essential  for  our  development.  It  is  the 
basis  of  each  technique:  decomposition,  partitioning,  and  basis  factorization. 

These  methods  have  been  used  extensively  in  large-scale  structured 
programming,  but  their  application  in  solving  stochastic  programs  has  been 
limited.  In  our  presentation,  we  demonstrate  how  these  techniques  can  be 
applied  to  the  stochastic  version  of  the  multi-stage  linear  program.  We  will 
then  show  that  our  algorithms  may  yield  substantial  savings  over  straight¬ 
forward,  “brute-force”  techniques  by  using  the  program  structure  effectively 
in  reducing  the  computational  cost. 

The  first  method  we  present  is  called  a  nested  decomposition  algorithm 


for  SDLP  or  NDSDLP,  because  it  decomposes  a  large  problem  into  succes¬ 
sively  finer  subproblems. 


The  basic  principles  of  optimisation  that  we  employ  are  outer  linearisa¬ 
tion  and  inner  linearization,  as  described  in  Geoffrion  [24].  Through  outer 
linearization  we  optimize  over  a  larger  convex  region  than  is  feasible  and 
then,  by  increasing  the  restrictions  obtained  from  the  subproblems,  approach 
optimality  within  the  true  feasible  region.  To  perform  inner  linearization  one 
optimizes  over  successive  subregions  of  the  feasible  region  and  approaches 
thereby  a  global  optimum. 

The  method  below  applies  outer  linearization  to  the  primal  problem 
of  SDLP.  It  can  also  be  viewed  as  applying  inner  linearization  to  the  dual 
problem.  We  then  have  a  master  subproblem  relationship  in  the  primal  as  in 
Benders’  method  [9]  or  in  the  dual  as  in  Dantzig- Wolfe  decomposition  [19]. 
We  also  follow  a  procedure  of  passing  between  periods  that  is  similar  to  the 
nested  decomposition  of  primal  inner  linearization  as  in  Glassey  [29]  and  Ho 
and  Manne  [32].  The  relationship  between  the  two  methods  is  well-known 
as  Kallio  and  Porteus  showed  in  [38]. 

We  begin  our  presentation  of  the  nested  decomposition  algorithm  in 
Section  2,  by  examining  some  properties  of  objective  functions  and  show¬ 
ing  the  basic  master-  subproblem  relationship  for  SDLP  under  the  outer 
linearization  scheme.  Section  3  presents  the  analogous  development  for  inner 
linearization,  or  Dantzig- Wolfe  decomposition.  We  follow  this  in  Section  4 
with  a  description  of  a  fundamental  problem  inherent  in  SDLP,  degeneracy. 
We  present  suggestions  for  its  resolution  and  some  of  its  special  difficulties  in 
the  stochastic  framework.  Lastly,  in  Section  5,  we  present  the  full  algorithm 
and  its  strategy  in  passing  through  the  scenarios. 
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i.  The  Master-Subproblem  Relationship 


The  decomposition  algorithm  we  discuss  below  relies  on  certain  fun¬ 
damental  properties  of  the  multi-stage  program  under  uncertainty.  These 
properties  concern  the  convexity  of  the  objective  function  and  the  repre¬ 
sentation  of  the  solution  set  as  a  convex  polyhedron.  Wets  first  observed 
these  attributes  in  [61]  and  [62].  These  are  reviewed  below. 

We  begin  by  formulating  the  equivalent  convex  program  to  SDLP.  We 
then  use  this  formulation  to  examine  how  we  can  find  the  set  of  linear 
constraints  that  constitute  the  solution  set.  This  procedure  involves  inducing 
feasibility  in  the  subproblems  and  finding  the  conditions  for  optimality.  The 
section  concludes  with  an  explanation  of  how  the  cutting  planes  for  these 
operations  are  constructed. 

Every  stochastic  dynamic  linear  program  (with  continuous  or  discrete 
distribution  of  the  random  variables)  is  equivalent  to  a  convex  program  with 
linear  constraints,  as  we  state  in  the  following  theorem  of  Wets: 

Theorem  1  The  stochastic  dynamic  linear  program  SDLP,  as  deGned  in 
Chapter  1,  is  equivalent  to  a  program  of  the  following  form: 

mincizi  +  Qi(xi) 

subject  to 
•Ai*i  =  » 

*i  6  D x,  ( ECP ) 

*i  >0, 

where  ,s  *  convex  function  and  is  a  convex  polyhedron. 


Proof.  See  Wets  [62].  | 
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The  Q  1(2:1)  of  Theorem  1  is  defined  as 

$i(*i)  =  E<2[Qi(*a,6)]  (1) 

where 

Q  1(221(2)  =  {min[c2X2+Q2(za)]|A2a:2  =  6  Z>2,  *2  >  0}.  (2) 

The  subsequent  Q*(z)  are  defined  iteratively.  We  can  solve  SDLP  by  repeated 
solutions  of  these  convex  programs,  but  the  objective  function  and  solution 
set  may  be  hard  to  find.  We  present  below  methods  for  finding  Qt{*t)  and 
Dt  without  explicitly  determining  the  functions. 

ECP  above  is  called  the  equivalent  convex  program  of  the  stochastic 
program.  We  note  that  this  theorem  holds  when  £  has  a  continuous  or  discrete 
distribution. 

Before  we  proceed  with  this  development,  we  further  note  that,  since 
every  period  of  SDLP  corresponds  to  an  identically  formulated  optimization 
problem,  we  can  form  a  subproblem  in  every  period,  t,  and  for  every  scenario, 

A  • 

j,  that  is  similar  to  ECP.  Let  j  be  the  immediate  ancestor  of  j,  and  q  = 
+  £},  we  have 

z\  =  mincta^  +  Qt(x3t) 

subject  to 

^  ^ 

Dizi  >  \l  ( ECP(t,j )) 

z{  >  0, 

where  we  have  written  the  convex  polyhedron,  D,  as  D{x\  >  and  Qt{z\) 
is  defined  as  in  (1). 


ECP(t,  j)  is  used  in  the  subsequent  analysis  to  deyelop  the  master  sub¬ 
problem  relationships  that  are  encountered  in  the  nested  decomposition  al¬ 
gorithm.  We  first  want  to  find  a  method  for  constructing  so  that  z\  will 
be  feasible  for  its  descendant  scenarios.  Let  J  be  a  descendant  of  j  in  period 
t  + 1  for  a  given  solution  x{'°  of  ECP (t,j).  We  have 

min^+1  =  et+  ix3t+1  +  Q(z^+1) 
subject  to 

At+ix?i+1  =  +  Btx{‘°, 

Di+1x{+1  >  X{+1.  (ECP(T  +  1,  J)) 

**+ 1  —  °» 

Now,  if  ECP(t  +1,7)  has  no  feasible  solution,  then  by  Farkas’s  lemma, 

r  t 


there  exists  a  vector 

such  that 

lOf+l-Ai+l  +2  l^t+i  ^ 

(3) 

(  2^4.1)  ^  0» 

(4) 

and 

l0H-i(ft-)-i  Btx\‘  )  +2  «rJ_|_1Xj_|_1  >  0. 

(5) 

So,  in  order  for  ECP(f  +  1,7)  to  have  a  feasible  solution  x\  must  be 
chosen  such  that 


— (l (T3t+1Bt)xi  >  +itrf+1#+1  +2^+iM+i-  (6) 

This  implies  the  following  lemma: 

Lemma  1.  For  every  descendant  scenario  of  j,  j,  in  period  t  +  1  and  for 
all  and  20*+1  satisfying  (3)  and  (4)  if  x\  is  feasible  in  SDLP,  then  it 

satisfies  the  inequality  (6). 
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Using  this  result,  we  can  add  (6)  as  an  additional  constraint  to  ECP(t,7). 
We  repeat  this  for  each  descendant  j  of  j.  We  solve  ECP(t,  j)  again  after 
all  of  these  cuts  have  been  added  and  proceed  downward  again  with  a  new 
value,  x3t'1.  In  this  manner,  we  iteratively  construct  the  constraint  set  for 
x{.  From  now  on,  for  clarity  of  exposition,  we  will  write  the  equations  of  the 
form  (6)as 

>  A+i((i+,h  (7) 

where  k  =  1, . . .  ,p. 

We  know  how  to  find  the  constraint  set  of  ECP(t,  j)  and  we  must  now 
find  a  method  of  constructing  the  convex  function,  Q(z$).  To  do  this,  we 
first  observe  that  ECP(t,y)  is  equivalent  to 

min  z{  —  ctz{  +  ${  (8) 

subject  to 

Atx{  = 

~(A+, Bt)z>  >  *  =  i . p, 

<5W)  >  9j, 

A  >  o. 

With  this  formulation,  we  consider  that  z3,0  is  again  a  solution  to  ECP(t,  j) 
and  that  ECP(t-(- 1,7)  is  feasible  for  all  J.  Next  let 

x{z3'°)J  =  ,T)\  *a(*t,0» D)  (9) 

be  the  optimal  dual  prices  vector  for  each  ECP(f-)-l,7),  given  z3t'°  a  solution 
to  ECP(t,  j),  and  define 

/.,W-°)  =  E<,+l[iri(4’0.j)-  ({?+,)] 


and 


to(xi’°)  =  Ee.+l(Z)  ■  (o^j^f+a)].  (11) 

*=i 

Now,  for  Jr(z{+1,  j)  optimal  in  ECP(*-f  1,  J)  for  any  z3t+l,  the  following 
inequality  must  hold 

+Btz\)+  £  *i(*i,J)<ri+2Cd,2)  > 

*=i 

*1  (*s‘°. DUi+l’,Btzi)  +  X)  *2(Xt>j)at+2((J+2>  (12) 

*=1 

and,  since  jt(x^  ,  J)  is  optimal,  we  have 

<3(4,  J)  =  + B‘4)  +  £  ■  «&(&.«)•  (i3) 

k=l 

Therefore,  by  taking  expectations, 

<?(*?)  >  ClW")  +  wW'°)  +  •  S.4-  (14) 

Letting  =  h(4’°)  +  Pzfaj'0)  in  (14),  we  have  the  following  lemma. 

Lemma  2.  If(zj,  ${)  is  a  feasible  solution  to  ECP(t,j)  written  as  in  (8),  then 

H  >  +  On (xt'°)Bt)x*. 

This  lemma  enables  us  to  form  additional  constraints,  as 

(k[  ■  Bt)z{  +  H  >pl  (15) 

for  successive  l,  where  k[  =  xi{z3‘l)  and  pl  =  p(z{'1).  These  cuts  are  then 
also  added  to  ECPft,/)  whenever  we  find  that  a  solution  (*f^0f’,)  is  such 
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that 


6{l  <  /  +  *[Btz{1.  (16) 

When  (15)  is  solved  for  all  k[,  and  pl,  then  we  have  achieved  master-subproblem 
optimality  between  ECP(i,  j)  and  its  descendants,  ECP(£,J). 

We  have  shown  how  the  master  problems  and  subproblems  can  be  con¬ 
structed  in  SDLP  by  using  the  fundamental  results  in  Lemma  1  and  Lemma  2. 
We  will  present  below  the  basic  algorithm  for  finding  master-suboptimality. 
This  algorithm  follows  closely  from  Benders  [9]  and  has  been  presented  in 
the  two-stage  case  by  Van  Slyke  and  Wets  [56].  It  is  an  outer  linearization 
scheme  because  the  feasible  regions  D\  and  convex  functions  are  succes¬ 
sively  approximated  by  the  inequalities  in  (7)  and  (15).  We  call  this  procedure 
OLSDLP,  for  outer  linearization  of  SDLP.  This  exposition  includes  Step  2', 
the  case  of  an  unbounded  solution  in  ECP(t,y).  The  justification  for  this 
procedure  can  be  found  in  Van  Slyke  and  Wets  [56]  for  a  deterministic  prob¬ 
lem.  We  omit  details  here,  because,  in  general,  we  will  use  the  algorithm 
with  upper  bounded  variables  and  no  unbounded  solution  will  be  possible. 

OLSDLP 

Step  1.  Solve  the  current  form  of  E CP (£,;),  using  Phase  I  and  Phase  II  of 
the  simplex  method: 

min  ctx\  -f- 

subject  to 
Atx{  — 

=  i ,  -,p, 

— (jt{  •  Bi)x{  +  e{  >  pl,t  =  1, . . .,q, 

A  >  o, 
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(17) 


where  we  set  p  =  q  =  0  initially  and  let  =  —  oo,  if  q  =  0.  If  (17)  is 
infeasible,  stop.  If  (17)  is  feasible  and  unbounded,  go  to  Step  2'.  If  (17)  is 
feasible  and  bounded,  go  to  Step  2. 


Step  2.  For  (z{‘x,  9%'*)  a  solution  of  (17),  solve  the  Phase  1  problem  of 
ECP(t  -f-  1,7): 

w  . 

minw*  =  ev  + 
subject  to 

At+lx3t+1  +lv  =  $+1  + 

^t+l*t+l  “t"  —  —  M+i 

■v 

xt+  i,v,u+,u~  >  0. 

(18) 

For  each  J,  such  that  wJ  >  0  in  (18),  use  the  resultant  multipliers  to 
build  a  cut  of  the  form  in  (7).  Add  these  cuts  to  (17)  and  increment  p.  If 
^>0  for  any  j,  go  to  Step  1;  otherwise,  go  to  Step  3. 

Step  2'.  From  (17),  we  obtain  an  unbounded  ray,  z{'1'  Xyj’1*,  for  X  >  0. 

Now,  solve  (18)  for  each  J,  but  replace  (J+a  +  Btz3,1  by  Bty\l* .  If  w?  >  0, 
for  any  j,  add  cuts  as  in  (7)  and  return  to  Step  1.  If  u3  =  0  for  all  j,  solve 
ECP(t  +  1,7)  for  all  J  with  the  same  replacement.  Let  Z*( t  +  1,  J)  be  the 
expected  value  of  the  objective  functions  and  compute  +  and  pl. 

Next,  solve  (18)  with  z3t'1  for  each  j.  If  ui3  >  0  for  any  j,  add  the 
feasibility  cuts  and  return  to  Step  1.  If  w3  =  0  for  all  J,  then  we  check  if 
ctv{'1  +  ^(t  +  1,7)  <  0.  If  so,  the  objective  function  is  unbounded,  stop.  If 
c*Vt’1  +  **(*  +  *»7)  >  0>  then  we  must  eliminate  yj’1  as  a  feasible  direction. 
We  do  this  by  using  the  icl  and  pl  found  above  in  forming  a  constraint  of  the 
form  (15)  and  adding  it  to  (17).  In  this  case,  we  return  to  Step  1. 
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Step  3.  Solve  ECP(H-1,  j)  for  each  j.  Compute  =  E^t+l 

•  ,  t  .  , 

x[,  and  px.  If  stop.  ECP(i,  ji)  is  solved.  Otherwise,  generate  a 

cut  of  the  form  (15)  add  it  to  (17),  and  return  to  Step  1. 

This  algorithm  terminates  in  a  finite  number  of  steps  as  we  state  in  the 
following  theorem. 


Theorem  2.  The  algorithm,  OLSDLP,  for  finding  the  solution  of  ECP(t,  j), 
results  in  either  an  infeasibility  criterion,  an  unbounded  solution,  or  an 
optimal  solution  in  a  finite  number  of  steps. 


Proof.  Every  iteration  of  the  algorithm  results  in  the  addition  of  a  constraint 
of  the  form  (7)  or  (15)  to  the  optimization  in  (17).  Since  there  are  a  finite 
number  of  bases  for  each  ECP(t  -f  1,  J),  the  number  of  these  constraints  is 
finite.  They  also  cannot  be  repeated  since  would  already  have  had  to 
satisfy  that  constraint.  Therefore,  the  algorithm  terminates  after  a  finite 
number  of  steps.  | 


The  finiteness  of  OLSDLP  can  also  be  maintained  if  we  delete  the  cuts 
in  (17)  that  are  slack  after  each  iteration.  This  is  true  because  the  objective 
function  in  (17)  is  monotonicaliy  decreasing  as  new  cuts  are  introduced.  We 
also  need  only  keep  at  most  m(t-f-l)+l  cuts  because  the  solution  of  ECP(f,y) 
has  at  most  m(t )  -f-  m(t-f"  l)-f- 1  basic  variables,  as  Murty  showed  in  [45].  We 
will  return  to  this  important  property  in  more  detail  in  Chapter  V  in  our 
discussion  of  the  local  basis  method  for  SDLP. 

We  have  seen  how  a  master  problem  at  period  t  in  SDLP  can  be  solved  by 
outer  linearization  using  subproblems  at  period  t-j~l-  As  we  stated  above,  this 
development  is  completely  analogous  to  applying  inner  linearization  to  the 
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dual.  In  the  next  section,  we  demonstrate  the  relationship  between  OLSDLP 
and  a  Dantzig- Wolfe  decomposition  form  of  inner  linearization. 

3.  The  Relationship  to  Dantzig- Wolfe  Decomposition 

Dantzig  and  Madansky  [18]  in  their  fundamental  paper  on  programming 
under  uncertainty  first  proposed  that  two-stage  stochastic  linear  programs 
could  be  solved  by  applying  Dantzig- Wolfe  decomposition  to  the  dual  of  the 
stochastic  linear  program.  In  the  context  of  our  development  here,  we  want 
to  apply  this  decomposition  to  the  linear  subproblem  of  SDLP  at  period  t 
and  scenario  j.  We  call  this  program  LP(£,  j).  (Note  that  in  this  case  we 
must  assume  that  £t+i  has  a  discrete  distribution  as  in  SDLP.)  The  problem 
we  address  is  then 


rninz^  =  etxi  +  + - 1-  p*et+iXt+i 

subject  to 

A,lJ  =  {j  +  B,_1xf_1 

Bt*(  +  —  Ct+i 

*  —  1»  •  •  •  1 

This  is  the  program  for  one  section  of  SDLP.  We  have  removed  the  constraints 
before  period  t  and  after  t  +  1. 

The  dual  of  LP(t,y)  can  be  written  as 

max  u(({  +  +plv1g£1  H - hp*w*^i 


subject  to 

uAt—p1vlBt  •••  pkvkBt  <  ct, 

<  Ct+i,  *  =  I,---.*- 


{DP(t,j) 
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Next,  let  {*■*}  =  {(jr^, . . . ,  xj^)},  where  /  =  1, . . . ,  q,  be  the  set  of  all  possible 
combinations  of  k  extreme  points  of  P  =  {x|jrAt+i  <  ct+i),  and  let  { am }, 
where  m  =  1, be  the  set  of  extreme  rays  of  P.  Define  also  T1  = 

P*jrkf,‘  and  lm,t  =  DP \j,t)  can  be  written 

in  Dantzig- Wolfe  decomposition  form  as: 

maxu(^'  +  Bi-ix3t_x)  -j-  XV  +  5Z  5Z  «'m,V/7m,‘ 


m=  1  t=l 


subject  to 


p  * 


-  E  *****  -EE  pm,ipiomBt  < 


m=  1  i=  1 


EX‘=l. 


/=rl 


X*  >  0;  pm,t  >  0,  for  all  t,  l,  and  m. 

(DWD(t,  j)) 

Now,  we  take  the  dual  of  DW(t,  j)  and  use  x3t  and  8  as  multipliers.  The 
result  is 

min  ctx\  -j-  8 
subject  to 

=  £ t  + 

{-ffmBt)x3t  >  »  =  1, . . . ,  *;  m  =  1, . . .  ,p;  (DWD(t  ~ 

( — rtBt)x{  +  >>/;  <  =  1 . (DW^,j)) 

z\  >  0. 

DWD(f ,  j)  has  the  same  form  as  (IT)  except  that  we  have  not  included 
the  extra  constraints  that  enter  into  the  subproblems  of  ECP(t  -j-  1 ,  J).  The 
feasibility  criteria  correspond  to  the  subproblems’  proposing  an  extreme  ray 
to  the  master  problem  in  DW(t,  j),  and  the  optimality  cuts  on  6  correspond 


to  extreme  point  proposals.  Optimality  in  the  outer  linearization  corresponds 
then  to  the  absence  of  better  proposals  from  the  subproblems  in  the  Dantzig- 
Wolfe  approach.  We  state  these  results  in  the  following  lemma: 

Lemma  3.  The  outer  linearisation  of  the  primal  problem  ECP{t,j)  in 
(17)  is  equivalent  to  solving  the  dual  of  Dantsig-Wolfe  inner  linearisation  as 
applied  to  the  dual  problem  of  ECP(t,j),  DP(t,j). 

Having  completed  the  analysis  of  this  basic  algorithm,  we  would  like  to 
show  how  it  is  implemented  in  solving  the  entire  program,  SDLP.  Further 
complications  enter  into  OLSDLP  because  of  possible  degeneracy  in  the 
subproblems.  In  the  next  section,  we  discuss  these  difficulties  and  how  they 
relate  to  stochastic  programs.  We  also  propose  ways  for  resolving  them. 

4.  The  Degeneracy  Problem 

One  weakness  of  decomposition  techniques  is  that  much  of  the  work  to 
optimize  subproblems  can  be  wasted,  because  the  final  inputs  from  the  sub 
to  the  master  differ  so  much  from  the  initial  ones.  This  can  lead  to  many 
iterations  from  master  to  subproblem  that  a  method  with  more  interaction 
between  the  problems  might  be  able  to  avoid.  The  next  two  methods  we 
present  have  a  more  unified  framework,  and  hence,  fit  this  description.  In 
this  section,  we  will  show  how  to  make  OLSDLP  more  responsive  in  the 
subproblems  to  changes  in  the  master. 

One  property  that  might  cause  unnecessary  iterations  in  decomposition 
schemes  is  the  fact  that  excess  columns  in  the  basis  of  the  master  problem 
(that  is, more  than  those  required  to  meet  the  original  set  of  constraints)  cause 
degeneracy  in  the  solution  of  the  subproblems.  Dantzig  and  Abrahamson 


[1]  first  observed  (and  then  proved)  this  property  in  their  experiments  with 
a  dual  nested  decomposition  algorithm  for  deterministic  multi-stage  linear 
programs.  They  also  noticed  that  repeated  sub-optimisation  changed  the 
basis  in  the  master  problem  very  slightly,  and  they  theorised  that  some 
efficiency  may  be  gained  by  allowing  the  subproblems  to  determine  some 
of  the  values  of  first  period  basic  variables.  This  is  possible  because  of  the 
subproblem  degeneracy. 

To  show  degeneracy  in  SDLP,  we  will  again  refer  to  ECP(t,  j)  as  written 
in  (17)  and  will  use  ECP(f  +1,7)  as 

-r  -r 

min^+i  =  ct+  i*t+1 


subject  to 


At+ iz{+1  =  +  Bti*, 

xj  >  0. 


(19) 

where  we  have  dropped  the  and  D\+x  from  ECP(t  +  1,7)  for  the  sake 
of  clarity.  Equation  (19)  is  the  form  of  ECP(t  +  1,7)  used  before  any  of  its 
subproblems  have  been  encountered.  The  degeneracy  result  is  included  in 
the  following  lemma. 


Lemma  4.  If  a  constraint  of  the  form  (7)  is  binding  at  the  optimal  solution 
x\‘* ,  to  (17),  then  every  feasible  primal  basic  solution  of  (19)  with  right-hand 

-r 

side,  +  Btx\’*,  is  degenerate. 


Proof.  For  the  binding  cut,  we  have 


0 


Now,  let  be  a  feasible  basis  for  (19).  By  applying  to  this,  we  find 

“H  &txi'  )  —  0-  (21) 

We  have  Vt+i^F+i  ^  °>  but  7^  0  for  a}^!  a  basis,  hence,  there 

exists  some  z¥+i,j  —  0,  proving  the  result.  | 

This  lemma  implies  that  a  degeneracy  will  occur  in  any  subproblem  that 
has  forced  a  tight  feasibility  cut  on  the  master  problem.  The  constraints  in 
0  which  enter  the  optimality  conditions  can  also  cause  degeneracies  in  the 
subproblems.  The  difficulty  with  these  degeneracies,  however,  is  that  they 
may  enter  in  any  subproblem,  and  we  may  not  be  able  to  determine  which 
one.  We  state  this  degeneracy  result  in  the  following  lemma. 

Lemma  5.  If  two  constraints  of  the  form  (12)  are  binding  at  the  optimal 
solution  *)  to  (17),  then  every  solution  of  (19)  for  all  j  which  satisfies 

T  . 

the  optimality  criterion,  *i+i  ^  includes  a  degenerate  solution  for  some 
1 • 

Proof.  Let  the  binding  constraints  be 


+  »{••  =  f\ 

(22) 

— +  *?>  =  p1. 

(23) 

Let  the  optimal  set  of  bases  for  {J1, ...  J*}  be  . . . ,  A^-*}.  Associated 

with  these  bases  are  prices  {jr1’*,...,**'*}.  These  prices  may  be  the  same 
as  those  for  one  of  (22)  or  (23).  Without  loss  of  generality  assume  they  are 
identical  with  the  prices  in  (23).  They  must  be  distinct  from  (22)  because, 


in  the  progression  of  the  algorithm!  (23)  must  have  been  violated,  when  (22) 
was  satisfied  and  before  (23)  was  added  to  the  constraint  set  of  (17).  Now, 
letting 


we  have 


and 


and/?1  =  53  P^tt+u 

»=i 

i=i 

k 

k 

=  E  + b,4 ••), 

(24) 

«=1 

t=l 

* 

k 

»=i 

=  Ep‘V,‘(«Vi  +  »**?'*)• 

i=  1 

(25). 

So,  bj  assumption 


t=i 


and  jt*>*  =  cBl,  and  <  cBi  for  jt1*1  feasible  for  all  t.  Therefore, 

we  obtain 

=  EWA  (») 


or 


^p*(c®*  —  —  0. 


(27) 


i=i 


Now,  pi  >  0,  x >  0,  and  eBi  —  <  0,  but  we  must  have 

eBi  x1,tAf^_l  for  jt*’*  unique;  therefore,  there  exists  some  ji  such  that 

“i . 

x3t+1  is  degenerate.  Hence,  the  result.  | 


The  degeneracy  we  have  shown  implies  that  the  subproblems  are  too 
restricted  to  alter  the  direction  of  the  solution  to  the  master  problem.  Dantsig 
and  Abrahamson  [1]  have  proposed  remedying  this  difficulty  in  the  multi¬ 
stage  deterministic  model  by  passing  columns  forward  from  the  master  to  the 
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subproblem  and  allowing  the  subproblems  to  determine  the  weights  of  these 
surplus  columns  in  the  master  problem.  We  present  below  an  application  of 
this  technique  to  SDLP  and  discuss  its  weaknesses  in  the  case  of  stochastic 
programs. 

We  assume  an  optimal  solution,  xf*,  of  (17)  is  partitioned  as 

*>■•  =  (*?,*;)  (28) 

where  xf  corresponds  to  a  square  nonsingular  submatrix  of  At,  A f ,  and  x* 
corresponds  to  A\.  We  can  write  xf  in  terms  of  x*  as 

I?  =  (ft® T\l\  +  ft+.*?+,)  -  (AfrlA^z-.  (29) 

So,  we  can  let  x*  vary  as  long  as  xf  >  0.  We  then  can  write  (19)  for 
each  j  as 

minct+1xj+1 
subject  to 

*.+!»?+ 1  =  {?+.+  Bf(Af )-'(({  +  ft-.  i?_.) 

*L-l  ^  0,  if  >  0, 


or,  letting 

*f  =  and  £+1  =  (t+l+B?(Af )-'((] ■ +ft_.i?_,), 


”T 

minct+ixf+1 

subject  to 

^r*f  =  £t+i> 

*t+i  —  0»  zf  ^  ^  0), 
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(31) 


where  we  hare  added  the  constraint  zf  >  0,  parenthetically,  because  we 
must  check  in  optimising  (31)  that  this  is  not  violated.  If  so,  we  would  fix 
the  variables  at  sero  and  proceed. 

This  problem  of  checking  for  feasibility  of  zf  enters  both  the  deter¬ 
ministic  and  stochastic  programs,  but,  in  the  stochastic  program,  any  of  the 
subproblems  of  type  (31)  may  determine  z*.  Since  z*  enters  the  program 
in  which  a  degeneracy  is  caused  (filling  the  degenerate  variable’s  position  in 
the  basis),  we  must  know  for  which  J  the  program  (31)  will  have  some  z* 
basic.  This  would  be  possible  if  all  the  z*  corresponded  to  tight  feasibility 
(type  (7))  cuts,  because  then  the  degeneracy  would  occur  in  the  correspond¬ 
ing  scenarios  which  generated  those  cuts.  For  the  optimality  cuts  (type  (12)), 
however,  degeneracy,  where  z*  is  basic,  can  be  in  any  scenario. 

To  apply  this  column  passing  technique,  in  general,  to  SDLP,  we  can 
formulate  the  following  alternative  form  of  (17).  It  includes  constraints  to 
keep  zf  feasible  as  well  as  the  remaining  additional  cuts. 

min  cf  zf  +  0$ 


subject  to 

-(Afr'Ajz;  >  -(Afr'iti + Bt-izLi)- 

* 

— (o|*k5*)z|  >  [ctk(£+1)],*  =  l,...,p, 

-(4Bl)x'  +  ei>p£,l=  l,...,q, 

>  0, 


(32) 

where  c*  =  ef(Af)~1All  -}-  e*,  and  the  other  quantities  under  tilde  are 
correspondingly  defined  to  reflect  the  substitution  of  (29)  for  zf . 

The  program  (32)  is  a  second  master  problem  that  we  can  use  to  deter¬ 
mine  the  optimal  values  of  z*  given  zf .  Our  proposal  then  is  to  follow 


OLSDLP  with  (32)  in  place  of  (17).  We  would  do  this  after  solving  the  sub¬ 
problems  of  the  form  (31),  which  have  generated  tight  feasibility  cuts  for  (17) 
and  in  which  we  know  degeneracies  must  occur. 

The  use  of  this  method  of  passing  some  columns  for  tight  feasibility  cuts 
and  then  determining  the  other  x\  by  using  (32)  in  OLSDLP  depends  on 
the  difficulty  of  solving  the  first  master  problem.  If  the  repeated  solution 
of  (17)  has  indicated  that  some  variables,  xf ,  are  persistent  in  the  basis, 
then  the  solution  of  (32)  could  obtain  the  optimal  values  without  involving  a 
reoptimisation  of  (17).  Furthermore,  if  the  number  of  surplus  variables,  x®  is 
sn  il,  (32)  may  become  significantly  easier  to  solve  than  (17).  This  alternative 
method  then  is  one  that  must  be  adopted  to  the  individual  problem  and 
its  requirements.  The  complications  of  creating  an  additional  optimization 
problem  in  (32)  may  outweigh  the  savings  in  solving  this  smaller  problem. 

5.  The  Complete  Methods 

We  are  prepared  now  to  present  the  algorithm  NDSDLP  for  the  entire 
stochastic  problem.  This  method  involves  repeated  use  of  the  algorithm 
OLSDLP,  and  proceeds  through  the  tree  of  possible  scenarios  in  SDLP  in  a 
forward  and  backward  manner.  Our  presentation  here  does  not  include  the 
resolution  of  the  degeneracy  problem  as  discussed  in  Section  4,  but  this  may 
be  added  as  a  modification  to  OLSDLP.  The  algorithm  follows. 

NDSDLP 

Step  0.  Set  up  a  problem  of  the  form  (17)  with  no  extra  constraints  for  each 
scenario  j  in  each  period  t  of  SDLP. 


Step  1.  Solve  (17)  for  period  T  (written  17-1).  Use  the  result  x\,  and  solve 
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(17)  for  each  node  in  period  2  ((17-2,;)  for  ;  =  1, . . . ,  *2).  If  any  subproblem 
is  infeasible,  add  a  feasibility  constraint  of  type  (7)  to  (17-1).  Solve  (17-1) 
again  and  proceed  to  period  2. 

After  each  subproblem  at  t  =  2  is  feasible,  proceed  to  solve  (17-3;)  for 
=  1, . . . ,  k3.  Again,  for  any  infeasibilities,  pass  back  cuts.  Continue  in  this 
manner  to  period  T,  until  there  is  a  feasible  solution  for  (17-t;)  for  all  t  and 
;.  If  an  infeasible  solution  to  (17-1)  ever  results,  then  stop-the  problem  is 
infeasible.  Else,  the  result  is  a  feasible  primal  solution  of  SDLP. 


Step  2.  Solve  OLSDLP  for  every  scenario  in  period  T  —  1.  This  implies 
master-suboptimality  for  the  last  period.  Set  T  —  T  —  2  and  go  to  Step  3. 

Step  3.  If  t  =  1,  go  to  Step  4.  Otherwise,  solve  OLSDLP  for  every  scenario 
in  period  t.  For  some  scenario  j  at  t,  this  may  involve  resolving  OLSDLP 
for  its  descendants  at  t  -j- 1  in  order  to  get  optimality  for  the  subproblems. 
We  say  a  subproblem  is  “solved",  in  terms  of  the  algorithm,  when  OLSDLP 
applied  to  it  ends  in  master-suboptimality. 

After  OLSDLP  has  been  solved  for  each  j  at  t,  set  t  =  t  —  1  and  repeat 
Step  3. 

Step  4.  Solve  OLSDLP  for  the  original  master  problem  at  period  1.  This 
program  may  again  involve  resolving  the  subproblems.  If  OLSDLP  at  1 
results  in  an  unbounded  or  infeasibility  termination,  then  stop-SDLP  is  ac¬ 
cordingly  unbounded  or  infeasible.  If  OLSDLP  ends  with  master-suboptimality, 
then  stop-the  current  solution  is  optimal  for  SDLP. 
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This  algorithm  follows  an  iterative  procedure  as  in  dynamic  program¬ 
ming.  We  pursue  this  relationship  more  closely  in  Chapter  VI.  NDSDLP  also 
terminates  in  a  finite  number  of  steps  as  we  state  in  the  following  theorem. 

Theorem  3.  The  method,  NDSDLP,  terminates  in  a  &nite  number  of  steps 
with  an  optima]  solution  to  SDLP  or  the  unbounded  or  infeasibility  conditions 
from  OLSDLP. 

Proof.  From  Theorem  2,  we  know  that  each  implementation  of  OLSDLP 
must  terminate  in  a  finite  number  of  steps.  Since  the  algorithm  proceeds 
backwards  after  each  period’s  scenarios  are  solved  by  OLSDLP,  the  terminal 
conditions  in  Step  4  must  be  met  in  a  finite  number  of  steps.  | 

Several  improvements  can  be  made  to  NDSDLP  to  aid  in  its  efficiency. 
We  have  already  mentioned  the  second  decomposition  possible  in  OLSDLP 
as  a  resolution  to  the  degeneracy  problem.  We  may  also  want  to  proceed  be¬ 
tween  the  periods  without  completely  satisfying  master-suboptimality.  This 
modification  may  help  efficiency,  but  it  must  be  done  carefully  in  order  to 
avoid  any  excessive  number  of  iterations  among  the  periods. 

Another  possibility  for  speeding  the  search  for  a  feasible  primal  solution 
in  Step  1  is  that  SDLP  may  have  inequality  constraints.  In  this  case,  the 
subproblem  at  period  t  is 

min  ctxt 
subject  to 

*T 

AtXt  >  (l  -f-  Bt—iXt—i, 

xt  >  0, 
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(33) 


for  all  j.  Hence,  to  find  feasibility  for  all  j,  we  need  only  solTe  the  Phase  I 
problem 

mineu 
subject  to 

AfZt  —  -fti  Iv  =  at  -j-  Rt—iXt—i 

Zt  >  0, 

(34) 

where  at  =  (a*.-  I  af  »•  =  maxT£^  .),  thus  at  >  &  for  all  j.  A  feasible 
solution  to  (33)  implies  that  each  subproblem  j  at  period  t  is  feasible  for 
Xt-l- 

By  solring  (33)  in  Step  2  of  OLSDLP,  instead  of  solving  the  Phase  I 
problem  (18)  of  each  j,  we  eliminate  many  unnecessary  optimizations.  This 
would  greatly  aid  the  efficiency  of  NDSDLP  for  SDLP’s  that  have  inequality 
constraints.  These  programs  are  quite  common  in  practice  and,  thus,  (33) 
should  prove  most  valuable. 

This  modification  and  our  presentation  of  NDSDLP  above  demonstrate 
some  of  the  possibilities  for  solving  SDLP  by  concentrating  on  the  optimiza¬ 
tion  of  smaller  subproblems.  In  Chapter  VH,  we  discuss  the  computational 
aspects  of  the  algorithm  more  carefully.  In  the  next  two  chapters,  we  present 
other  methods  that  rely  upon  subproblem  optimization,  but  maintain  closer 
ties  between  the  sub-  and  master  problem. 
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Chapter  IV 
A  Piecewise  Strategy 


1.  Introduction 

The  decomposition  algorithm,  NDSDLP,  in  the  last  chapter  broke  the 
stochastic  program  SDLP  into  a  sequence  of  master-subproblem  relations. 
A  drawback  to  NDSDLP  could  come  from  wasted  effort  in  optimizing  sub¬ 
problems  without  affecting  the  master  problems'  inputs.  The  approach  in 
this  chapter,  the  piecewise  strategy  for  SDLP,  or  PCSDLP,  also  separates  the 
program  into  master  and  subproblems,  but  it  allows  for  only  one  optimization 
of  the  subproblems  without  the  master’s  involvement.  This  method,  which 
maintains  some  ties  among  the  separate  scenarios,  forms  a  bridge  between 
the  decomposition  approach  in  Chapter  III  and  the  local  basis  factorization 
of  Chapter  V. 

The  method  is  called  “piecewise’  because  it  relies  on  the  piecewise 
linear  property  of  the  objective  function.  Piecewise  methods  in  general 
(see  Geoffrion  [27])  follow  an  optimizing  trajectory  across  the  regions  of  the 
feasible  set.  For  a  convex  function,  an  optimization  is  performed  on  each 
region  that  leads  either  to  a  boundary  or  interior  solution.  If  the  solution 
is  interior,  then  that  point  is  optimal.  If  a  boundary  point  is  optimal,  one 
optimizes  on  the  adjacent  region  and  repeats  the  process.  (See  Figure  1.)  If 
no  direction  in  an  adjacent  region  is  improving,  then  the  current  point  again 
is  optimal. 

The  piecewise  strategy  has  been  applied  to  large-scale  linear  program¬ 
ming  through  a  method,  called  “partitioning”,  for  which,  J.  B.  Rosen  [53] 
has  been  most  responsible.  Our  use  of  the  strategy  in  this  chapter  will  be  to 
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exploit  the  repetitions  of  the  blocks  of  SDLP  and  to  form  an  algorithm  that 
can  adapt  to  different  scenarios  and  combine  them  adequately.  In  Section  2, 
we  present  the  basic  master-subproblem  algorithm  for  this  method  and  show 
where  efficiencies  can  be  made  in  its  implementation.  Section  3  then  states 
the  strategy  for  the  full  program  and  presents  the  difficulties  that  may  occur 
with  PCSDLP. 

2.  The  Master-Subproblem  Relationship 

The  method  proceeds  by  performing  two-period  optimizations  for  suc¬ 
cessively  larger  problems.  In  this  section,  we  present  the  two-period  problem, 
in  which,  the  first,  parent,  scenario  forms  the  master  problem  for  its  direct 
descendants,  the  subproblems.  Without  loss  of  generality,  we  assume  that 
this  optimization  takes  place  between  the  first  and  second  periods. 

The  first  period  problem  we  solve  is 


'/Igure  1.  The  piecewise  path  begins  at  z0  and  leads  to  zs. 


min  f°(*i)  —  ei*i 

subject  to  AiXi=  6if  (1) 

*i>  0. 

If  (1)  is  infeasible,  the  program  is  infeasible,  and  we  stop.  If  (1)  is 
unbounded,  then  we  follow  Step  2*  of  NDSDLP  to  remove  this  case.  If  we 
succeed,  we  return  with  the  cuts  that  eliminate  the  unbounded  raj  and  resolve 
(1).  Now,  for  x°,  an  optimal  solution  to  (1),  we  want  to  find  Q(z?)  as  ECP 
of  Chapter  m,  where 

Q{x°)  =  E$#[  min  62*2  s.t.  A2x2  =  f2  +  B\x\,  x2  >  0].  (2) 

Here,  if  there  exists  f£  such  that  there  is  no  feasible  solution  in  the 
jth  scenario,  then  we  form  a  cut  as  in  (3.7)  and  add  it  to  (1)  as  part  of 
its  constraint  matrix.  We  continue  until  each  subproblem  is  feasible.  Next, 
associated  with  each  there  exists  an  optimal  basis,  A$ ,  for  the  problem 
in  <?(z°).  We  write  Q(x%)  as 

*» 

«(*?)  =  E  (3) 

The  function  Q(x j)  from  (3)  is  linear  for  all  x2  feasible  for  Thus, 
fci 

«(*i)  =  E  0  « 

j=i 

for  all  x  such  that 

*2 = >  0.  (*<>) 


where  we  have  assumed  a  single  lower  bound  for  all  x2  at  0.  In  general, 
and,  for  most  practical  purposes,  both  lower  and  upper  bounds  should  be 

67 


included,  augmenting  the  constraints  in  (4a).  In  implementations,  this  is 
especially  important  for  artificial  variables,  where  lower  and  upper  bounds 
coincide. 

Now,  we  can  write  f°(*i)  in  (1)  as 

C°(*i)  =  («.  +  £  (5) 

3=1  j=l 


for  all  x\  such  that  +  (A%*)~ 1Bxxi  >  0  for  all  j  =  1, . . . ,  *2. 

We  note ,  for  optimal  multipliers  ic*  in  subproblem  j,  that  jrJ  =  c%s 
which  makes  computations  in  (5)  somewhat  easier.  We  also  note  that  each 
need  not  be  unique,  so  we  may  use  one  basis  several  times  without  check¬ 
ing  (4a)  for  each  j.  We  will  return  to  this  idea  below,  but,  for  the  current 
exposition,  we  will  use  constraints  (4a)  explicitly.  Thus,  for  x\  restricted  as 
in  (4a),  we  look  for  an  improving  direction  along  the  current  linear  section 
of  Q(*i)  by  solving 


min  f l(il)  =[a  +  EjLx  P’eatd?')  lBx]xx 
subject  to 


Aixx 

z 


=  h,  (60) 

>^,j  =  l,...,*a,(66) 

>0, 

(«) 


where  p3'  —  (A®')-1^- 

From  (6),  we  obtain  an  optimal  solution,  x\,  and  optimal  objective  value, 
fl(zj).  This  is  the  initial  point  of  our  feasible  path.  We  next  look  at  the  set 
of  constraints  in  (6b)  which  are  active.  We  define 


T  =  {(•,/) :  IWf'r'BiK'.*)'*!  =  ft, 


(T) 


where  p{  is  the  «th  component  of  the  vector  fP .  T  includes  the  rows  and 
scenarios  that  generated  a  binding  cut.  It  corresponds  to  the  set  of  degeneracies 
in  the  subproblems  of  NDSDLP  and  is  directly  related  to  the  surplus  columns 
described  in  the  next  chapter.  We  next  order  the  pairs  (t,  j)  £  7  lexicographi¬ 
cally  and  write  7  =  {ri,  fa, ,  rp}. 

Now,  if  T  =  0,  then  z°  is  a  solution  of  (6),  so,  by  our  convexity  result 
in  Lemma  2.1,  *(z$)  =  eiz?  ~f  <J(z$)  <  ctz i  +  Q{x i)  for  all  feasible  Z\. 
Hence,  z®  is  an  optimal  solution  of  this  two-stage  SDLP. 

For  7  7^  0,  we  consider  subproblem  j  for  ri  =  (*,  j).  Here,  we  have 

+ (aZt'b^  m 


and 

*?'(•-)  +  =  0,  (9) 

where  x£  is  the  non-basic  partition  of  the  variables,  z2.  We  now  want  to 
force  zf'(t)  out  of  the  basis,  so  we  can  follow  a  new  path  in  the  adjacent 
feasible  region.  To  maintain  optimality,  we  find  the  entering  variable  as  in 
the  Dual  Simplex  Method,  (see  Dantzig  [16]) 


min 


(10) 


-rtf,.  ,  -rtf,.  ' 

— Za  (s,  a)  -X,  (*,;)> o  ^  (t,  j) 

where  and  are  the  representations  of  c$  and  relative  to  the  basis, 

We  can  now  pivot  out  xf*(i)  and  replace  it  with  z^(s).  In  doing  this, 
we  keep  and  call  the  new  basis  ,  an  auxiliary  basis,  where 


(11) 


for  17(1,  s),  the  elementary  matrix  corresponding  to  the  pivot  of  z^(s)  into 
the  basis  in  position  t.  This  new  basis  is  used  to  restrict  Xi  in  (4a). 

We  now  can  formulate  our  auxiliary  problem,  for  /t  =  1, 


min  {>■•(*,)  =[t,  +  £*!, 
subject  to 


A 1X1=  b  i, 


(12) 


(^)_1SiXi>  pf’\S=l,...,ka, 


zx>  0. 


In  (12),  we  hare  changed  the  objective  function  from  (6)  and  some  of  the 
constraints  (6b).  The  variables  z\  still  form  a  feasible  basis,  but  they  may 
no  longer  be  optimal. 

We  proceed  to  price  out  the  cost  row  in  (12)  with  the  new  parameters 
and  to  check  optimality.  If  z\  is  still  optimal,  then  we  drop  and  return 
to  T  for  r2,  again  find  an  entering  variable,  and  form  the  auxiliary  problem 
(12).  As  long  as  we  cannot  improve  on  z\  for  rlt  we  try  tj+i.  If  we  find 
/  +  1  >  p,  then  no  direction  can  improve  on  z{,  hence,  it  is  optimal. 

If  we  find  that,  for  any  rj,  z\  is  not  optimal  in  (12),  then  we  optimise 
(12)  and  obtain  z\.  We  set  fa(zf)  =  f1,#(z^),  replace  A%*  by  form  a 

new  set,  T,  of  the  tight  constraints,  and  proceed  again  with  rx  in  search  of 
am  improving  direction. 

We  follow  the  above  procedure  to  obtain  a  sequence  of  decreasing  objec¬ 
tive  values,  f4(z()  >  ^(z?)  >  •••  >  ?M(zf),  until  we  cannot  improve  our 
current  solution.  The  algorithm  we  have  described  has  the  following  steps  : 


PCSDLP(2) 

Step  0.  Find  an  optimal  bounded  feasible  solution  of  (1),  or  terminate  if 
infeasible  or  if  Step2l  of  NDSDLP  implies  unboundedness. 
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Step  O'.  Find  an  optimal  feasible  basis,  ,  for  each  subproblem  /  in  Q(x\) 
by  applying  feasibility  cuts  (3.7)  to  (1). 

Step  1.  Form  the  program  (6)  using  the  set  of  bases,  B  =  Solve  (6) 

and  obtain  and  x\.  Form  7.  If  T  =  0,  stop,  x°  is  optimal.  If  T  5^  0, 
set  l  =  1,  p  =  1,  and,  go  to  Step  2. 

Step  2.  If  l  >  p,  stop,  x£  is  optimal.  For  r<  =  (t,y),  find  the  entering  variable 
a  in  scenario  j  by  (10).  Form  the  auxiliary  basis,  and  the  program 

(12)  for  fM,*(xi).  Using  as  a  starting  solution,  solve  (12)  and  obtain  x\.  If 
fM'*(xJ)  =  fM’*(*i),  set  l  =  /+1  and  return  to  Step  2.  If  fM,*(xJ)  <  ?M,*(xf), 
go  to  Step  3. 

Step  3.  Update  T  and  B.  Set  =  xj,  l  =  1,  p  =  p  + 1,  and  go  to  Step 
2. 

The  following  theorem  states  the  finiteness  of  PCSDLP. 

Theorem  1.  The  method  described  above,  PCSDLP,  terminates  in  a  Unite 
number  of  steps  with  an  infeasible,  unbounded,  or  optimal  solution  to  the 
two-stage  (T  —  2)  form  of  the  program,  SDLP. 

Proof.  From  Chapter  HI,  we  know  that  Steps  0  and  O'  must  terminate  in  a 
finite  number  of  steps.  After  Step  2,  the  solution  to  (6)  and  (12)  must  be 
feasible  in  SDLP  because  primal  feasibility  of  the  last  solution  is  maintained. 
It  is  bounded  because  (6)  and  (12)  are  more  restrictive  than  (1). 

T  is  finite  since  the  number  of  constraints  (12b)  is  finite  and  each  im¬ 
proving  solution  corresponds  to  a  new  set  of  bases,  B.  Since  there  are  a  finite 
number  of  possible  basis  set  combinations,  the  algorithm  must  terminate.  | 


3.  A  Method  for  Redueed  Basil  Storage  Requirements 

Efficiency  and  storage  requirements  in  the  solution  of  (12)  can  be  significantly 
improved,  if  we  do  not  include  redundant  constraints  that  occur  with  dupli¬ 
cated  bases,  .  To  do  this,  begin  by  checking  in  Step  O'  for  a  repetition  of 
A^.  We  start  with  /  =  1  and  increase  l,  letting  each  distinct  new  basis  be 
Aj *.  We  obtain  B  =  {A%l\l  =  1, . . . ,  q).  We  also  define 

us) 

j(0 

where  J(l)  =  {  all  scenarios  j  with  optimal  bases,  Af*}. 

Now,  when  we  construct  the  constraints  (12b),  we  define 


/>'-*=  max  jeJ{i)P3,*» 


(14) 


and,  for  each  component  s’,  we  store  }'(/,*)  for  every  l,  where  ^'*(t')  > 
for  all  j  €  J{1).  Thus,  (12)  becomes 


min  f**(*i)  =[ci  +  ]C?=iP(0c?'’*(A?‘’*)  lBi]zx 
subject  to 


A\X  j=  &i, 

(A?'*)"1^  1=1, 

Xj>  0. 


(15a) 

•7,(156) 

(15) 


Now,  in  Step  2,  write  the  elements  of  T  as  i\  =  ($, }'(/)).  Each  time  a 
new  auxiliary  basis,  A^'*,  is  investigated,  if  6  B,  then  we  adjust  p(l) 
and  possibly  but  do  not  change  the  coefficient  matrix.  If  $  B, 
then  we  must  add  another  set  of  constraints  to  (15).  At  Step  3,  we  update 
B,  the  associated  probabilities,  p(t),  and  update  T  using  the  pair,  (s,/(Q). 
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This  modification  can  significantly  reduce  the  number  of  constraints 
since  a  large  number  of  scenarios  may  have  the  same  basis.  We  could  again  use 
the  Garstka-Rutenberg  procedure  (see  Chapter  I)  to  find  the  probability  that 
each  basis  is  optimal  in  Step  O'  without  solving  the  individual  problems.  This 
effect  combined  with  the  smaller  size  of  (15)  can  lead  to  greater  computational 
efficiencies. 

Another  efficiency  can  be  gained  from  using  a  method  similar  to  the 
column  passing  technique  of  Chapter  m.  During  the  algorithm,  we  may 
observe  that  one  set  of  variables,  {zf  },  remains  in  the  optimal  basic  set, 
{zf},  while  the  other  variables  are  chosen  from  a  set,  {zf }.  If  the  columns 
of  the  set  {zf }  have  full  rank,  we  can  take  a  square  non-singular  submatrix, 
Af ,  from  (15)  and  find 


*?  =  (A?)"1#,  -  (A?)-*A? xf .  (16) 

We  then  eliminate  zf  from  (15)  and  obtain 

min  f**(zf )  =  c^zf 

subject  to 

-Kz  f  >  -(Af)-1*!, 

1 . q, 

*?  >  0, 

(17) 

where  tilde  indicates  that  2J,  A*,  and  p1'*  are  defined  relative 

to  zf  through  substitution  of  (16)  into  the  program  (12).  The  definitions  are 
completely  analogous  to  those  in  (3.32). 

The  optimisation  procedure  can  then  continue  with  the  reduced  problem 
(17)  to  find  the  optimal  values  zf,  given  zf .  Using  the  result  of  (17)  in 
(15)  would  then  determine  optimality.  When  decisions  can  be  narrowed  to 
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choices  among  a  few  variables,  this  modification  may  again  prove  effective  in 
improving  efficiency. 

The  algorithm,  PCSDLP(2),  is  stated  for  two-period  optimisations.  It 
does  not  require  the  subproblems  to  be  reoptimized  after  Step  O'.  By  main¬ 
taining  primal  feasibility,  it  is  always  on  a  feasible  path  to  the  solution  and 
may  eliminate  the  problems  of  wandering  among  multiple  suboptimal  points. 

In  the  next  section,  we  present  the  implementation  of  PCSDLP(2)  for  general 
multi-stage  programs. 

4.  The  Complete  Solution  Strategy 

The  PCSDLP  method  follows  a  procedure  very  similar  to  NDSDLP  in  its 
passing  through  the  scenarios  from  period  to  period.  In  fact,  both  of  these 
methods  can  be  seen  as  local  approximations  of  a  dynamic  programming 
scheme,  which  we  present  in  greater  detail  in  Chapter  6.  PCSDLP  even 
begins  by  finding  a  feasible  primal  solution  through  NDSDLP,  but  PCSDLP 
never  allows  for  primal  infeasibility  or  non-optimality  in  a  subproblem  after 
a  single  optimisation. 

First,  we  set  up  subproblems  for  each  node  as  in  Step  0  of  NDSDLP. 
Next,  we  find  a  feasible  primal  solution  by  Step  1  of  NDSDLP,  passing 
feasibility  cuts  as  we  proceed  through  the  periods  from  1  to  T.  After  finding 
this  feasible  solution,  we  start  by  applying  PCSDLP(2)  to  the  master-subproblem 
relations  at  period  T  —  1.  Primal  feasibility  is  then  maintained  throughout 
the  optimisation. 

For  each  scenario  j  in  period  T  —  1,  we  solve 
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min 

subject  to 

At-ix3t_1  =  4— 1  +  Bt— a*r_2» 

-r  -r  _  _ 

— Bt— \x3t_x  +Ar*^=  4»  j  =  1, . . *r, 

4-i  ^  °.4>  °»  ?= 

by  using  PCSDLP(2). 

For  a  solution  to  (18),  define 

T  T 

P3  =  {• :  Xj.(»)  is  basic  in  subproblem  j) 


(18) 


and 

-r  __ 

T5  =  {* :  *  is  basic  in  row  I  and  (l,j)  E  T). 

Now,  with  the  solution  from  PCSDLP(2),  we  want  to  find  the  optimal 
basis  for  the  full  problem  (18).  This  larger  matrix  will  form  the  basis  of  a 
subproblem  for  period  T  —  2.  We  first  include  the  set  of  basic  variables  in 
the  master  problem,  The  basic  variables  from  the  subproblems  will 

be  chosen  as 

X’r  =  {4(.'):ie/3?r|7?}  (19) 

This  definition  eliminates  the  degenerate  variables  from  the  basic  set.  Since 
the  elements  of  X3T  are  the  only  non-zero  variables  in  an  optimal  feasible 
solution  to  (18),  if  (18)  is  not  degenerate,  then  the  union  of  and  X3T 

must  form  a  basic  set  of  variables  in  (18). 

If  (18)  is  degenerate,  then  we  must  check  whether  the  columns  cor¬ 
responding  to  {zbtLx)  X3T  span  the  solution  space.  If  not,  we  add 
columns  from  those  corresponding  to 

Tr  =  (4(.)  :  i  £  4 Hi7)  (20) 
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If  a  column  of  xJT(t)  £  Xy  is  independent  of  the  columns  in  the  present  basis, 
then  it  is  added  until  full  rank  is  a  achieved.  We  know  that  we  must  obtain 
a  basis  since  the  union  of  all  columns  for  { x X3T,  and  Xy  spans  the 
space. 

We  thus  obtain  a  basis  for  all  in  T  —  1.  We  call  this  basis,  Dp_lt  and 
we  further  define 


Dt-i  = 


dr— 1  = 


(  At—  i 
— Bt—i  At 

\—Bt—  i  At  ) 

(cr—  \,P1ct,--,P*tct), 


and 


(21) 

(22) 

(23) 

(24) 


4-i  =  (&-..  4.  •••.<rr)r. 

4-i  =  (4-i>  4‘ . *tt) 

Now,  the  problem  for  PCSDLP(2)  is 

min  f£_2  =cr-a4-a  +  EjLT  P3* r-il^-i  (25) 

subject  to 

At-2Z*t_2=  Z3t-2  +  Bt—\z3t—V  (25o) 
—Bt-2Z3t_2  +  Dt- i=  i>3r_ v  j—  1,  — ,  — 1,(256) 

XT—  2  —  Vt—  1  —  3 '  ~  1»  •  •  • 

where  Bt—2  i*  the  matrix  Sr— 2  augmented  by  zeroes  to  correspond  with 
&T-V 

To  begin  PCSDLP(2)  for  (25),  we  substitute  constraints  of  the  form  in 
(6b)  for  (25b)  and  enter  Step  1  in  PCSDLP(2).  We  never  reoptimize  the 
subproblems,  but  look  for  feasibility  maintaining  pivots  in  both  periods  T — 1 


and  T.  In  general,  after  solving  the  two-stage  problem  for  each  scenario 
j  in  period  t,  we  would  again  find  the  bases,  Df * ,  and  construct  a  two- 

-T 

-r  Jj. 

stage  problem  for  t  —  1  as  in  (25)  by  combining  x\  with  (j/^, . . . , yf£+l). 
PCSDLP(2)  would  begin  by  optimizing 

min  =[c-i  + 

subject  to 

■A-t — — 1=  ff—x  +  Bi — 2xt—2> 
(Dfi)~1Bt—iX3t_1  >  p3 ,  j  = 

*Li>  o. 

The  dual  pivoting  operations  could  then  be  performed  in  any  of  periods  t 
through  T. 

PCSDLP  continues  by  combining  master  problems  with  subproblems 
and  following  these  iterations  back  to  period  1.  This  process  uses  the  basis 
structure  depicted  in  Figure  2b.  of  Chapter  I,  in  which,  we  view  each  scenario 
as  starting  a  new  problem.  The  steps  we  have  described  for  PCSDLP  follow. 

PCSDLP 

Step  1.  Follow  Step  0  and  Step  1  of  NDSDLP  to  obtain  a  feasible  solution 
t,o  SDLP.  Set  t  =  T  —  1  and  set  up  a  program  of  the  form  (12)  for  each 
scenario  j  in  T  —  1.  Set  j  =  1. 

Step  2.  Follow  Steps  1,  2,  and  3  of  PCSDLP(2)  for  the  problem  at  node  (j,  t). 
If  j  =  ktl  go  to  Step  3.  If  j  <  kt)  set  j  =  j  +  1  and  return  to  Step  2. 

Step  3.  If  t  =  1,  stop,  zf  i*  optimal.  If  /  >  1,  combine  the  master  and 
subproblems  of  each  scenario  j  at  period  t  and  form  programs  as  in  (26)  to 
initiate  PCSDLP(2).  Set  t  =  t  —  1  and  go  to  Step  3. 

The  finite  termination  of  this  method  is  guaranteed  by  the  finiteness  of 


PCSDLP(2)  and  our  passing  back  one  period  in  each  encounter  with  Step  3. 
We  state  this  as  a  Corollary  to  Theorem  1. 

Corollary.  PCSDLP  terminates  in  a  finite  number  of  steps  with  an  infeasible 
or  unbounded  solution  from  t be  procedures  of  NDSDLP  or  with  an  optimal 
bounded  feasible  solution  to  SDLP. 

PCSDLP’s  greatest  potential  improvement  over  NDSDLP  is,  as  we  have 
emphasized,  its  maintenance  of  primal  feasibility  and  subproblem  optimality. 
This  advantage  over  the  possible  suboptimizations  and  infeasibilities  in  NDSDLP 
must  be  discounted,  however,  by  the  growth  of  the  bases,  Df ,  in  the  sub¬ 
problems.  Their  larger  size  may  lead  to  a  greater  number  of  computations 

in  performing  pricing  and  the  minimum  ratio  test  in  (1).  We  can,  however, 
gain  efficiency  with  a  compact  factorization  of  Df.  The  following  chapter 
describes  such  a  technique  and  its  application  to  the  full  problem,  SDLP. 
This  local  basis  method  could  then  be  used  in  conjunction  with  PCSDLP  to 
gain  still  greater  efficiency. 
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CHAPTER  V 


A  Local  Baaii  Simplex  Method 


1.  Introduction 

The  two  methods  presented  above,  NDSDLP  and  PCSDLP,  require  the 
optimization  of  subproblems  essentially  independent  of  the  main  problem. 
This  practice  can  be  costly,  leading  to  a  great  number  of  iterations.  We 
describe  below  an  adaptation  of  the  Simplex  Method  for  linear  programs  with 
stochastic  structure.  This  method  reduces  the  complications  of  stochastic 
linear  programs  by  taking  advantage  of  some  fundamental  properties  of  the 
basis.  We  call  this  approach  a  local  basis  simplex  method,  LBSMPX  , 
because  it  relies  on  the  near  square  block  triangularity  of  the  bases  for  these 
problems. 

Block  triangular  linear  programs,  in  general,  have  the  form: 

min  c  1X1  -J-  c 2X2  -{-  •  •  •  -f-  ct%t 
subject  to 

(1) 

AtyXi  H - h  AttXf —  bt,  t  —  1,  •  • .,  T, 

xt>0,t=  1, . . . ,  T. 

where  xt  €  R^fit  €  Rm*, ct  €  Rn*,  and  the  matrices  AtJ  are  dimensioned 
accordingly. 

The  detached  coefficient  matrix  is  then: 

fMx  \ 

Mi  Mi 


A  = 


Tl  At2  •  •  •  AtT 
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(2) 


rp  tp 

and  has  dimension, m  X  »,  where  m  =  ]£t=i  »nt  and  n  =  Ylt=.i  ni- 

Dantzig  [15]  introduced  the  idea  of  using  an  artificial  basis  for  these 
programs  in  which  square  ( m *  X  mt)  basis  blocks  appears!  along  the  diagonal. 
The  true  basis  could  be  derived  from  the  artificial  basis  by  a  set  of  side 
conditions.  He  observed  that,  because  of  the  persistence  property  in  dynamic 
linear  programs,  the  additional  computations  would  be  few. 

This  concept  of  basis  factorization  and  the  use  of  pseudo-  basic  variables, 
as  in  Beale [6],  have  been  applied  in  a  variety  of  examples.  Recent  implemen¬ 
tations  are  found  in  Kallio  and  Porteus[38],Perold  and  Dantzig[48],Fourer[22], 
and  Propoi  and  Krivonozhko[52].  This  last  approach  is  close  to  the  develop¬ 
ment  here,  using  local  bases  as  a  factorization  for  multi-stage  linear  programs. 
A  thorough  and  unified  presentation  of  the  relationships  among  basis  fac¬ 
torization,  partitioning,  and  decomposition  can  also  be  found  in  Winkler  [67]. 

2.  The  Structure  of  the  SDLP  Basis 

The  program, SDLP,  is  another  member  of  the  class  of  block  triangular 
linear  programs  as  the  following  shows. 

Lemma  1.  The  program,  SDLP,  has  the  structure  of  a  block  triangular 
linear  program,  as  in  (1). 

Proof.  Define 


(Ax 


Att  — 


Ax 


\ 


\  At) 


(3) 


where  Ax  is  repeated  kx  times  to  correspond  with  each  scenario  in  period  t. 


Define  also 


-Bt-x 


-Bt-x) 


(4) 


where  kt—x  repetitions  of  the  matrix,  —Bf—x,  correspond  to  the  possible 
outcomes  of  period  t  —  1. 

Then,  we  define  At_iit  as 


The  other  Ay  matrices  are  void,  so  SDLP  has  the  desired  structure.  | 

The  SDLP  has  other  advantages  that  it  shares  with  general  multi¬ 
stage  linear  programs.  In  these  programs,  the  number  of  additional  columns 
required  in  finding  the  true  basis  from  the  artificial  basis  is  limited.  This 
bound  means  that  storage  requirements  should  not  grow  excessively  with  the 
problem  size.  This  well-known  property  is  stated  in  the  following  lemma: 

Lemma  2.  For  a  dynamic  linear  program  (  a  block  triangular  linear  program 
where  A,*  =  0  for  all  a  >  t  -j-  1,  ie.  a  staircase  structure),  the  number  of 
surplus  and  deficient  columns  in  each  block  of  the  basis,  Bu,  is  bounded  by: 
0)  0  <  lx  <  m2, 

(ii)  — m*  <  lt  <  =  2, . . . ,  T  —  1, 

(Hi)  b  >  -mf. 

where  lt  -(-  is  the  number  of  basic  columns  in  the  basis  from  period  t,and 
the  basis  has  the  form: 


Bn 

S21  B22 

0  B&2  S33 


/ 


A 


B  = 


(5) 


l  0  0  0  Btt—i  Btt* 


Proof,  (ia).  In  order  for  B  to  be  nonsingular  Bn  must  have  full  rank. 
Therefore,  l\  >  0. 

(ib)  and  (iia).  lt  <  ntt+i.  Assume  lt  >  rrit+i,  then  B  has  greater  than 
mt  -J-  mt+i  independent  columns  in  period  t,  but  the  row  rank  of  this  block 
is  less  than  or  equal  to  mt  -f-  a  contradiction.  Hence,  It  <  mt+i. 

(iib)  and  (iii).  /*  >  — m<.  Assume  lt  <  —mt,  then  there  are  greater 
than  ms  independent  columns  in  periods  1, 2, . . . ,  t  —  1,  which  again 

contradicts  row  rank,  so  lt  >  — m*.| 

The  direct  application  of  this  lemma  to  the  stochastic  linear  program, 
SDLP,  would  imply  that  for  any  period,  t,  there  are  at  most  kt  -mt+i  surplus 
columns  in  the  basis.  Because  of  the  highly  structured  nature  of  SDLP, 
however,  a  much  tighter  bound  can  be  found.  We  show  this  below  in  the 
following  lemma.  Murty  [45]  first  observed  this  property  for  the  two-stage 
case  of  programming  under  uncertainty  with  continuous  distributions  of  the 
random  variables  in  1968. 

Lemma  3.  For  SDLP  and  for  each  scenario  j  in  period  t 

(i)  0<h<  m2, 

(ii)  —mt  <  l{  <  mt+i,j  =  =  2,...,T—  1, 

(iii)  —mt  <  l{;j  =  1, . . . ,  kT. 

where  l\  is  the  number  of  surplus  columns  in  the  basis  for  period  t  and 
scenario  j,  and  mt  is  the  number  of  rows  for  a  single  scenario  in  period  t. 
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Proof,  (ia)  This  follows  from  Lemma  2. 

(ib)  and  (iia).  Assume  that  1$  >  mt+1.  This  implies  that  there  are 
greater  than  m*  -f  independent  columns)  in  the  matrix  : 


(X\ 


(6) 


\-Bj 


but,  again,  Dt  has  row  rank  mt  +  ntt+i,  a  contradiction.  Therefore,  < 
r*t+ 1* 

(iib)  and  (iii).  As  before,  the  number  of  independent  columns  cannot 
exceed  the  row  rank  of  the  submatrix  including  that  scenario,  and  the  result 
follows.  | 

From  this  result,  we  proceed  to  find  an  efficient  implementation  of  the 
Simplex  Method  to  the  problem,  SDLP,  in  which,  a  limited  number  of 
additional  computations  are  used  to  perform  the  simplex  routines  of  finding 
the  true  basis  representation  of  a  column,  the  cost  row,  and  prices. 

3.  Finding  the  true  basis  representation  of  a  column 

To  find  the  column  representation,  first  define  a  square  block  triangular 
artificial  basis  for  SDLP  as 
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Vi 

Vi 

n 

Vi 

Ul 

Vh 

(7) 

Vi 

• 

V  VTT-\  UTT' 

where  U3t  consists  of  the  rows  of  A3  for  mt  independent  columns  in  each 
period  t,  and  Vt  includes  the  corresponding  rows  in  — Bt.  For  a  constraint 
matrix,  A,  we  also  define 

a  =  {aj  E  A  :  ay  is  a  column  in  the  artificial  basis}. 

We  further  define 

t  =  {dj  E  A  :  Oj  is  a  column  in  the  true  basis  }. 

The  complements  of  a  and  r  are  defined  as  V  aad  T. 

For  every  U\,  there  also  is  a  partition 


\P’t  T\) 


where  P\  and  T\  are  square  and  non-singular,  the  columns  in  P{  are  pseudo- 
basic,  chosen  from  a  (")  T,  and  T\  contains  columns  from  the  true  basis.  We 
define  here 

=  {*  :  Uj  E  a  Qf  where  a j  is  the  basic  column  for  row  Ik}. 

To  replace  these  columns  in  the  basis  at  period  t,  there  must  be  surplus 
basic  variables  from  period  t  —  1.  We  write  these  columns  in  terms  of  the 
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artificial  basis  as 


where  Pf-i  consists  of  rows  from  AJ—1  and  is  m*_  i  X  s{_1  ( the  number  of 
surplus  columns  in  period  t  —  1  under  scenario  j  ).  The  rows  corresponding 
to  — are  partitioned  between  an  a3t_l  X  a\_^  matrix,  and  R\_v 
a  (*t  •  mt  —  «^_i)  X  *t_i  matrix. 

The  coefficient  matrix, A,  is  then 


which,  by  row  and  column  permutation,  is 


The  following  allows  us  to  use  this  partitioning  for  computations. 

Lemma  4.  Relative  to  the  artificial  basis,  U,  the  true  basis,  T,  of  SDLP  has 
columns  partitioned  as  in  (11),  where  each  QJt  is  square  and  nonsingular. 

Proof.  By  the  definition  of  the  surplus  columns  in  (9)  and  suitable 
permutations,  we  arrived  at  (11). 

Next,  define  Q  as  the  matrix,  relative  to  the  artificial  basis,  of  surplus 
columns  with  rows  in  V>>  that  is 
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which  is  clearly  nonsingular  because  U  1  and  T  are  both  nonsingular. 

Now,  index  Q\  for  all  j  and  t  along  the  diagonal  as  Qx, . . . ,  Qk.  Let  Qp 
be  the  first  singular  matrix  in  the  list.  Its  row  rank  is  rp  <  ap,  the  number 
of  columns.  In  order  for  Q  to  be  nonsingular,  however,  the  row  rank  of  the 
submatrix  of  the  remaining  columns,  Fp+i,  is 

k 

Fp+i  >  ai‘ 


This  violates  the  column  rank,  therefore,  each  Q3t  is  nonsingular. | 


Given  that  the  are  nonsingular,  we  can  proceed  to  eliminate  the 
pseudo-basic  columns  from  the  basis.  This  procedure  involves  premultiplying 
U~XA  by  a  product  of  matrices,  F*! I71.  We  first  define 


PiQT1 

RiQT1 


I 


/ 


(12) 


and  observe  that 


F1U~1A  = 


/  Pi 


?2 


We  next  define 


*1  = 


and,  in  general, 


F\  = 


Q  r 


— i 


QkT—l 

r — 1 


—  i 


-«i9r‘ 


'i  -vim-1  1 

(«a)-‘ 

V  V,(Ql)  1 

'/  -Piwi)-1  \ 


(13) 


(14) 


wir1 

V  — H1,  (<?{)->  ij 

where  Ht  and  pj  are  the  partitions  of  the  surplus  columns,  relative  to  the 
previous  FJg’ s.  By  repeated  multiplication,  the  matrix  relative  to  the  true 
basis  is  found: 

<1 

Fjfi  •  •  F\U~XA  —  [  F^  -Fi  FkTT-1i---F1U~1A 

I 


(15) 


In  order  to  facilitate  finding  the  matrices,  F3g,  we  note  that  the  growth 
of  nonseroes  in  every  period  and  scenario  is  limited.  The  following  lemma 
states  this. 

Lemma  5.  After  multiplication  by  Flt~ l,  F\~2,  ...,F  J, . . .,  F\,  the  additional 
noniero  blocks  in  the  surplus  columns  of  scenario  l  in  period  t  occur  only  in 
rows,  for  which,  I’s  ancestor  scenarios  are  basic. 
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Proof.  We  proceed  by  induction  on  p,  where  we  order  all  the  blocks  of  surplus 
columns  with  p  =  1, . . .,  k.  For  p  =  1,  the  result  is  true  trivially  since  no 
new  nonzero  blocks  are  added. 

Assume  this  is  true  for  all  p  <  n.  We  look  at  the  surplus  columns  in 
n  + 1,  and  assume,  without  loss  of  generality,  that  this  is  scenario  l  in  period 
t.  We  then  define 


Next  observe  that  every  Fp  in  periods  1  —  2  has  only  identities  in 

columns  corresponding  to  P\,  Q\,  and  R\,  so  that  only  the  Fit_l  need  be 
considered.  If  the  scenarios  are  not  ancestors  of  l  at  period  t,  then  again  by 
the  hypothesis,  they  have  no  block  entries  that  correspond  to  Rlt,  Qlt,  or  P\, 
and  so  do  not  alter  S\. 

Therefore,  we  have 


and  the  only  additional  blocks  occur  where  —  is  nonzero. 

By  the  hypothesis,  the  only  nonzero  blocks  in  are  in  its  ancestor  rows 
from  previous  periods,  so  the  only  additional  nonzero  blocks  in  Slt  will  occur 
in  these  rows  and  the  rows  where  columns  are  surplus  in  t  —  1,  scenario  l. 
This  completes  the  induction.  | 


This  lemma  proves  valuable  in  computing  the  columns  relative  to  the 
true  basis  .  Multiplications  and  storage  are  limited  since  only  sections  of 
each  surplus  block  need  to  be  considered.  We  next  wish  to  compute  the 
value  of  a  column,  Zi,  where  *,•  6  T(t,  l),  the  non  -basic  columns  in  period  t 
and  under  scenario  l. 

We  first  observe  that 

Fn--F1U~1zi  =  FltFTt_1? 7=27*,  (18) 

where  7  is  the  first  ancestor  scenario  of  l,  since  for  every  F?  such  that  j  £ 
{t,t—  1}  orp^i 

FjZ7  =  27  (19) 

by  the  definition  of  F?. 

Next  partition  the  components  of  the  vector  Zi  as 

S?  =  fri :  j  e  m) 


6?  =  (ty  :  j  6  1>(t  +  1)) 

K  =  ■  3  e  m) 


K  —  (2y»  :  j  G  +  !))•  (20) 


Then,  z\  can  be  written  as 


1 

1 

►-* 

£ 

1 

»-» 

>-* 

»J  = 

(Q\- 1)"1 

•*?  + 

«{r* 

-HUiWUi)-1 

l  J 

-site!)-' 

t 

K 

\  j 

(21) 


With  this  representation,  we  form  the  following  procedure,  called  LBFTRN 
(for  local  basis  forward  transformation)  to  find  Z* . 

LBFTRN 

Step  0.  Identify  an  incoming  column  z%.  Find 


Z.  =  ((^n1af;(l/;+ir1(Vi((C/ir1a<  +  6,))  =  (*.;5.)as  defined  in  (20). 

(22) 

Step  1.  Identify  3?  and  Jf.  Find 


(23) 


Step  2.  For  all  j  G  U^«=i  iH*)  (*M  preceding  pseudo-basic  columns)  find 


(24) 


for  j  e  V»(t  ~  1), 


2jf=W-P[(j,  *)F, 


(25) 


for  /  G  V>(t), 


i 


" j * 


for  j  G  ULi  *(*)> 


Zj,  =  *)w  ~  R\U>  *)v> 


and  for  j  G  i>{t  —  1), 


Step  S.  For  j  G  0, 


rit- — Rte*)®- 


(26) 

(27) 

(28) 
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(29) 


and  for  j  6  i>{t  +  1,0 


(30) 


Zji  =  Z3<  +  Ki- 

LBFTRN  results  in  the  updated  vector, I*.  Therefore,  from  (21)  and  our 
previous  results,  the  following  proposition  holds. 

Proposition  1.  The  representation  of  an  incoming  variable,  z%,  relative 
to  the  current  basis  in  SDLP  can  be  found  as  IS*  from  the  above  procedures 
in  LBFTRN. 

This  routine  allows  for  quick  computations  of  an  incoming  column  when 
there  are  few  surplus  columns.  It  requires  current  information  of  which  rows 
belong  to  ip(t)  and  1),  the  pseudo-basic  columns  in  period  t  and  period 
t  -f- 1,  the  representations  of  P\_ 1(  P\,  R\_ lf  and  the  inverses  of  Q\—i 

and  Qlt.  This  storage  requirement  would  be  small,  relative  to  the  storage  of 
the  entire  matrix,  if  few  surplus  columns  are  present. 

4.  Finding  the  Dual  Prices  and  Pricing 

The  backwards  transformation  involved  in  computing  the  dual  prices  for 
SDLP  can  also  be  done  more  efficiently  with  the  use  of  a  square  block  trian¬ 
gular  artificial  basis.  A  savings  here  could  be  substantial  because  the  pricing 
operation  aften  requires  a  majority  of  the  effort  in  linear  programming  codes 
(viz.,  Fourer  [22]  ).  The  method  presented  below  requires  only  the  present 
and  previous  scenario  blocks’  surplus  column  updates  for  computation  of  the 
reduced  costs. 

We  assume  again  that  we  are  in  period  t  and  at  scenario  l.  The  dual 
prices  relative  to  the  artificial  basis  are  computed  first,  by  transformations 
from  the  last  period  to  the  beginning.  We  define 

? rip  =  c ^.(£7^)  xfor  all  /  =  1, . . .,  kr- 
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(31) 


1 


and  define  xj  by  the  general  recursion: 


for  all  l  =  1, . . .,  kt,  where  l(t  -}-  1)  denotes  all  descendants  of  l(t)  at  t  +  1 
and  is  the  probability  of  each  descendant. 

Now,  to  find  the  prices  relative  to  the  true  basis,  we  look  for  p  such  that 

0  =  Cj  —  p  •  aj  for  all  j  E  r,  a,  E  A.  (33) 

We  have  from  jt  defined  in  (32)  and  (33)  that 

0  =  Cj  —  w  •  a j  for  all  j  E  r  P|  a, 

or,  since  jt  =  cU~l, 

0  —  Cj  e(U  1oy) 

=  Cj  —  Cj  •  Ij, 

where  Ij  is  a  unit  column  with  identity  in  row  j. 

Thus,  for  c  partitioned  as 

c9  =  (c j  :  j  E  ^P)  a),  and 


(34) 

(35) 


cR  =  (cj.jerf |  a),  (36) 

we  have  from  (35)  and  the  definition  of  Q  and  R  in  (9), 

h  =  c:  ~  5Z  c9i '  pi  -  cQi '  Qi  ~  c*  ’  Ri  *>r  all  y  E  r  P|ff.  (37) 

•yy 

We  seek  next  a  such  that 

0  =  Cj-  £(e?  +  ai)Pj  -  (c«*  +  Oj)Qj  -  cRRj  for  all  j  E  rpjff.  (38) 
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We  do  this  iteratively  by  finding  first 


0 1  =  1 


(39) 


and,  in  general, 

=  [tfj  -  for  all/  =  1, . . .,  kt.  (40) 


We  define  p  as 


if  j€l> 
if  j  6 


(41) 


Hence,  (33)  holds  by  the  construction  of  a  in  (39)  and  (40)  and  the 
preservation  of  eR  =  0. 

The  computation  of  p,  using  a  as  defined  in  (39),  allows  pricing  in  an 
individual  period  t  to  involve  only  the  inverse  of  Qlt.  If  pricing  proceeds  from 
period  t  to  t  -f  1,  then,  for  a  constant  artificial  basis,  the  previous  prices 
need  not  be  recomputed.  This  strategy  may  result  in  substantial  savings  by 
eliminating  unnecessary  multiplications  for  the  other  periods.  It  also  saves 
on  storage,  since  extraneous  surplus  block  inverses  can  be  ignored. 

The  local  basis  simplex  procedures  for  finding  dual  prices  are  then: 


LBBTRN 


Step  0.  For  period  t  and  scenario  l,  cr  =  (cy(f,  l) :  j  6  rf)5)- 
Step  1.  Find  6j(t,  l)  as  defined  in  (37)  for  all  j  6  r{t}  l)  f)  5. 
Step  2.  Compute  a\  as  in  (40). 

Step  3.  Form  p  as  in  (41). 

This  procedure  is  then  followed  by  LBPRCE  which  finds 


Ij  =  cj-p •  a, 
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(42) 


for  all  j  £  T{t,  l).  The  incoming  variable  in  this  strategy  is  then  chosen  as 
the  variable  with  the  most  negative  reduced  cost  in  that  period  and  scenario. 
We  define  this  as 

?,(t,  l)  =  mini€r(t)0Cj  -  p  ■  <ij.  (43) 

If  Zi(t,t)>0,  then  the  algorithm  would  proceed  to  the  next  scenario  in 
period  t,  or,  if  l  =  kt,  the  first  scenario  in  period  t  +  1  would  be  considered. 

5.  Updating  the  Pseudo-Bases  and  Surplus  Bloeks 

After  the  incoming  column  is  chosen  in  (43)  a  leaving  column  is  selected 
by  the  minimum  ratio  criterion  of  the  standard  simplex  method.  This  pro¬ 
cedure,  called  CHUZR,  finds 

=  min,,.>0(Sl),  (44) 

where  is  the  representation  of  the  incoming  column  found  by  LBFTRN 
and  is  the  current  value  of  the  right  hand  side  in  the  ith  row. 

The  computations  in  CHUZR  can  also  be  reduced  because  nonzero  entries 
in  zj  are  restricted  to  certain  blocks  as  we  discussed  in  Section  3.  For  up¬ 
dating  the  basis,  s  can  be  pseudo-  basic  or  non-basic  in  period  t  (ie.,s(t,  l )  £ 
a  Pj  T(t,  l)  or  s(t,  l)  £  Q  y(t,  /)  ),  and  s  can  enter  the  true  artificial  basis 
in  period  t  or  become  surplus  basic  in  period  t  -f-  1  (that  is,  r  £  r(t  — 
1, l)  n*,r(f  +  l,/)n <*>  -f  1, l)  f| ff;°r  r(t  +  1, 0  D <*)■  We  discuss  these 
cases  below  : 

Ul.  r£r(f— 1,/). 

In  this  case  »  replaces  a  surplus  column  from  the  previous  period.  To 
update,  we  remove  a  column  and  row  from  Q{t  —  1,  l)  and  update  the 
corresponding  list  of  pseudo-vectors  in  period  t.  We  also  update  U(t,  [)~1 
by  replacing  the  pseudo-basic  column,  that  corresponded  to  r  with  s. 
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This  is  done  by  adding  an  elementary  matrix  to  the  eta  file  of  U(t,l)~ 1  as 
in  standard  simplex  codes.  (Note  that  if  a  £  a,  the  artificial  basis  remains 
unchanged.) 

U2.  r  £r(t, Of) <*• 

Here,  s  replaces  a  currently  true  basic  column  in  the  artificial  basis.  If 
8  £  a,  then  we  must  replace  the  row  in  Q(t,  l )  for  which  a  is  basic  with  the 
row  for  which  r  was  basic.  The  artificial  basis  is  unchanged.  If  a  £  at,  then 
we  need  only  update  the  artificial  basis,  U(t,l),  as  in  Ul. 

U3.r  6  r(t  -f- 

In  this  case,  s  replaces  a  surplus  column  which  is  basic  in  the  next  period. 
We  maintain  the  same  artificial  basis  and  update  the  surplus  block, Q(t  + 
1,  /).  To  do  this,  the  column  occupied  by  r  in  Q(t  +  1,  /)  is  replaced  by  the 
corresponding  row  entries  of  a.  This  is  performed  easily  by  premultiplying 
by  an  elementary  matrix  for  this  pivot. 

U4.  rer(t+l,0fl«- 

In  this  case,  s  replaces  a  basic  column  in  the  next  period.  This  involves 
adding  a  column  and  row  to  Q(t  -f- 1,/).  The  new  row  has  entries  from  each 
of  the  surplus  basic  columns.  The  list  of  pseudo-basic  variables  in  the  next 
period  must  also  be  updated  with  the  addition  of  the  row  for  which  r  was 
basic.  Again,  the  artificial  basis  remains  unchanged. 

The  updating  procedure  can  be  confined  to  the  current  local  artificial 
bases  and  current  and  following  pseudo-bases.  This  property  enables  us  to 
store  only  the  present  and  following  bases  for  updating  purposes. 

0.  The  Algorithm 

The  previous  sections  have  presented  the  basic  routines  for  a  local  basis 
simplex  method.  We  discuss  below  the  method’s  basic  strategy  and  im- 
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plementation.  The  algorithm  is  similar  to  the  nested  decomposition  approach 
of  Chapter  III  in  that  it  alternately  follows  a  forward  and  backward  procedure 
through  the  periods.  This  strategy  is  also  similar  to  dynamic  programming 
since  the  algorithm  proceeds  from  one  local  optimum  to  the  next.  The  basic 
method,  call  LBSMPX,  follows. 

LBSMFX 


Step  0.  Find  an  artificial  basis.  Do  this  by  solving  first 

min  t\Z\ 
s.  t.  Axxx=  £, 

*i>0, 


and  proceeding  to  solve 


min  ctz\ 

s.t.  Atzlt—  (i  +  ft-tijf*"1), 

*i>0, 


(45) 


(46) 


for  all  (t,  /),  t  =  l,...,T,l  =  l,...,fc«.  If  any  of  these  programs  has  no 
feasible  solution,  use  the  last  infeasible  solution  basis  as  U\.  Otherwise,  we 
store  each  locally  optimal  basis  as  lflt  and  form  the  artificial  basis  from  these 
matrices. 

Step  0.  Set  t  =  T  -  1,1  =  l.NDRCTN  =  ‘ BACK'.CFLAG  = 
4 NO' ,M  ST  AT  =  'YE  S’, IT  NO  =  0. 

Step  1.  Invert  the  current  basis  and  find  all  basic  variables  values.  This 
procedure,  called  INVERT,  is  equivalent  to  performing  LBFTRN  on  the  right 
hand  side.  If  all  variables  have  feasible  values,  MSTAT  =  'YES’.  Else, 
MSTAT  =  ‘NO’. 


Step  2. Form  a  Phase  I  objective  row,  if  MSTAT  =  'NO’.  Else,  use  the 
objective  row,ct.  This  routine  is  FORMC. 
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Step  3.  Call  LBBTRN  for  (t,l)  to  find  the  current  prices.  Set  ITNO  = 
J7WO  +  1. 

Step  4.  Call  PRICE  for  (t,l).  If  7,  >0,  check  basic  variable  values.  If 
there  exists  z(t,l/,)  infeasible,  set  MSTAT  =  ' NO ’,  and 

(a)  If  t  =  T,l  =  kr.NDRCTN  -=  'FORE', and  CFLAG  =  ‘NO’, 
if  MSTAT  =  ‘ YES' ,  the  solution  is  optimal, stop.  If  MSTAT  —  ‘ NO * 
and  the  objective  value,  z(l,  1)  >  0,then  the  solution  is  infeasible,  stop.  If 
1(1,  l)  =  o,  then,  set  MSTAT  —  ‘ YES’ ,  and  go  to  Step  2. 

(b)  If  t  =  1,/  =  1  ,NDRCTN  =  ‘ BACK ’,  and  CFLAG  —  'NO',  then 
if  MSTAT  =  ‘YES',  the  solution  is  optimal,  stop.  Else,  if  the  objective 
z(l,  1)  >  0,  then  there  is  no  feasible  solution,  stop.  If  z(l,  1)  =  0,  then  set 
MSTAT  =  ‘YES'  and  go  to  Step  2. 

(c)  If  t  —  T,l  —  kT,NDRCTN  =  ‘ FORE and  CFLAG  =  ‘y£S’,set 
*  =  T  —  1,1  =  1,NDRCTN  =  ‘BACK'.CFLAG  =  lNO\ set  MSTAT  = 
y2?S.  If  ITNO  <  MAXIT[the  maximum  number  of  iterations  between 
reinversions), go  to  Step  2,  else  go  to  Step  1. 

(d)  If  t  =  1,/  =  1,IVI>RCTN  =  ‘BACK’,  and  CFLAG  =  TES’, 
then  set  t  =  2,1  =  l.NDRCTN  =  ‘BACK’.CFLAG  =  ‘NO’, MSTAT  = 
‘YES’.  If  ITNO  <  MAX  IT,  go  to  Step  2.  Else,  go  to  Step  1. 

(e)  If  t  <  1  and  /  <  kt,  then  set  l  =  l  -f-  1.  If  ITNO  <  MAXIT,  go 
to  Step  2, else  go  to  Step  1. 

(f)  If  1  <  f  <  T,l  =  kt.NDRCTN  =  ‘BACK’,  set  t  =  t—  1,1  =  1,  go 
to  Step  2. 

(g)  If  1  <  t  <  T,l  =  ktlNDRCTN  =  ‘FORE’,  set  t  =  t  +  1,/  =  1, 
go  to  Step  2. 

If  7,  <  0,  set  CFLAG  =  ‘YES’. 

Step  5.  Call  LBFTRN. 


Step  6.  Call  CHUZR.  If  no  r  is  found  (ie.,  all  ?*,<0),  then  the  solutnn 
is  unbounded,  stop. 

Step  7.  Call  UPDATE  to  perform  Ul,  U2,  U3,  or  U4,  depending  on  r, 
go  to  Step  2. 

The  preceding  method  leads  to  an  optimal,  unbounded,  or  infeasible 
solution  as  the  following  theorem  states. 

Theorem  1.  Assuming  nondegeneracy  (or  suitable  resolutions  by  choice  of 
outgoing  variable  in  CHUZR  (see  Dantzig  [17])),  the  method,  LBSMPX, 
terminates  in  a  finite  number  of  steps  with  an  optimal  solution,  an  unbounded 
feasible  solution  or  an  infeasibility  criterion. 

Proof.  Each  pass  through  Step  7  results  in  a  decreased  objective  value  under 
nondegeneracy.  Since  there  are  a  finite  number  of  bases,  the  algorithm  can 
pass  Step  7  at  most  a  finite  number  of  times.  The  method  would,  therefore, 
terminate  as  in  the  simplex  method  as  long  as  it  takes  at  most  a  finite  number 
of  steps  between  passes  of  Step  7. 

The  algorithm  does  not  return  to  Step  7  for  each  iteration  that  Step 
4  results  in  (c),  (d),  (e),  (f),  or  (g).  However,  in  these  cases,  the  algorithm 
progresses  to  the  next  scenario  or  period  in  either  the  forward  or  backward 
phase.  Therefore,  conditions  (a)  or  (b)  in  Step  4  must  be  met  after  each 
phase. 

When  (a)  or  (b)  are  encountered  in  Step  4,  the  algorithm  proceeds  in  the 
opposite  direction,  if  Step  7  was  ever  encountered  in  the  last  pass  through 
every  scenario  ( the  case,  CFLAG  =  ‘YES’).  If  no  improvement  was  made, 
that  is,  if  Step  7  was  not  passed  in  the  last  phase,  the  algorithm  terminates. 
Therefore,  LBSMPX  always  terminates  in  a  finite  number  of  steps.  | 

As  in  the  Simplex  Method,  variations  in  LBSMPX  may  be  used.  The 


99 


most  negative  pricing  strategy  in  PRICE,  for  instance,  can  be  altered.  We 
can  also  choose  to  proceed  to  the  next  period  with  a  criterion  other  than 
?,(£,  /)>0.  ?«(*,  l)>  —  e  could  be  used  with  c  decreasing  in  size  as  the  optimal 
solution  is  approached. 

LBSMPX  uses  the  standard  simplex  techniques  by  efficiently  partition¬ 
ing  the  basis.  As  in  the  nested  decomposition  approach,  NDSDLP,  and  the 
piecewise  method,  PCSDLP,  LBSMPX  concentrates  on  different  sections  of 
the  basis  one  at  a  time.  It  differs  with  the  previous  methods,  however,  by 
continuously  reflecting  changes  in  the  entire  problem.  The  advantage  of 
this  property  is  that  the  global  solution  reflects  the  local  optimization  more 
quickly.  The  disadvantages,  however,  are  that  the  method  requires  more 
computational  effort  in  maintaining  these  instantaneous  changes.  These  rela¬ 
tive  computational  requirements  are,  again,  discussed  more  closely  in  Chapter 
VII. 
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Chapter  VI 
The  Relationship  to 
Dynamic  Programming 


1.  Introduction 

In  Chapter  I,  the  dynamic  programming  formulation  for  the  general 
stochastic  dynamic  linear  program,  SDLP,  was  presented.  In  principle,  SDLP 
could  be  solved  exactly  by  this  method.  Since  the  random  vectors,  £t,  were 
allowed  to  have  continuous  distributions,  the  program  at  stage  t  was  an 
infinite  dimensional  optimization  problem.  To  avoid  this  complication,  we 
approximated  the  distribution  with  discrete  valued  random  vectors,  fj,  and 
formulated  a  linear  program,  the  solution  of  which  we  discussed  above. 

We  considered  linear  programming  techniques,  but  we  could  have  used 
dynamic  programming  to  solve  this  discrete  distribution  problem.  In  a  produc¬ 
tion  example,  Beale,  Forrest,  and  Taylor  [7]  estimated  that  four  state  vari¬ 
ables  in  a  time  period  could  be  handled.  For  larger  problems,  they  concluded, 
an  approximation  must  be  used. 

We  chose  the  linear  programming  formulation  because  of  the  histori¬ 
cally  good  performance  of  the  Simplex  Method  and  simplex-based  algo¬ 
rithms.  Using  these  algorithms,  large  problems  with  over  a  thousand  rows 
and  columns  can  be  solved  easily  (see  White  [66]  for  a  discussion  of  computa¬ 
tional  experience  with  large-scale  linear  optimization).  Dynamic  program¬ 
ming  techniques  using  the  standard  computational  procedure  are  generally 
limited  to  problems  with  six  state  variables  and  six  decision  variables  (see 
Larson  and  Casti  [41]).  Advanced  techniques  can,  however,  be  implemented 
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for  larger  problems  in  order  to  approximate  and  then  refine  the  state  and 
policy  space.  These  methods  prove  very  valuable  for  general  transition  equa¬ 
tions  and  objective  functions  (see  Larson  [39]),  but  linear  programming  methods 
are  most  commonly  used  for  problems  with  linear  constraints  and  objectives. 
They  do  not  require  quantization  of  the  state  space  or  knowledge  of  the 
range  of  state  variable  values  along  the  optimal  path.  Linear  programming  is, 
therefore,  often  more  efficient  than  dynamic  programming-type  algorithms. 

Our  methods  in  Chapter  IE  and  IV  actually  combine  these  two  tech¬ 
niques,  although  we  have  presented  them  as  linear  optimization  strategies. 
In  the  following  sections,  a  standard  dynamic  programming  formulation  of 
the  stochastic  linear  program  will  be  presented  and  its  advantages  will  be 
discussed.  We  will  then  show  that  the  piecewise  method,  PCSDLP,  and  the 
nested  decomposition  approach,  NDSDLP,  are  simply  different  versions  of 
general  dynamic  programming. 

2.  The  Quantized  State  Space  Approach 

A  standard  simplification  in  dynamic  programming  is  quantization.  SDLP 
can  be  formulated  in  this  manner  by  first  creating  a  discrete  approximation 
of  the  state  space.  This  will  enable  us  to  form  a  recursion  at  every  stage  t 
and  for  every  state  yt.  We  first  define 


yt  =  Btxt.  (1) 

and  quantize  this  vector  as  {yj,  y*;  •  • . ,  yj*}.  The  random  right-hand  sides 
will  again  have  discrete  values,  {fj,  £* , . . . ,  £**}. 

Now,  the  following  backward  algorithm,  BDP,  can  be  used. 
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BDP 


Step  0.  For  all  y*T,  find 

zt(v't-i)  =  E£T[*r(Vr— 1>  £r)]»  (2) 

where 

«(y‘r-ii^r)  =  min  cTxT 
subject  to 

(3) 

Atzt=  +  £r> 
xt>  0. 

Set  t=T.  Go  to  Step  1. 

Step  1.  Set  t  —  t  —  1.  For  all  find 

(4) 

where  6)  =  min  ft),*?)},  where 

z*((vi- v &)> Vt)  =  min  c<*t  +  *t+i(v?) 

subject  to 

■Atzt=  ^  +  W~i.  (5) 

—Btxt=  yf, 
xt>  0. 

Step  2.  If  t=2,  go  to  Step  3.  Else,  go  to  Step  1. 

Step  3.  Find  zx  =  zx(0, $i)  =  min  yyZi((0,&i),p{),  where  «i((0,6i),|^)  is  as 
defined  in  (5)  for  yo  =  0,  (i  =  b\.  Stop. 

Proposition.  The  algorithm,  BDP,  converges  in  a  finite  number  of  steps  to 
an  optimal  solution,  for  Zt(vj-i),  at  ever^  sta8e  t  ot  the  PTogrzm 

SDLP,  under  the  following  assumptions  : 
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I.  The  support  of  the  random  rector,  (t,  is  B*  =  {fj,  £*,...,  £*}. 

II.  The  values  of  the  inventories  from  period  t  are  B%Zt,  where  B\Z\  £ 

Yt  =  {*«>&•• -,1ft}* 

Proof.  The  result  follows  directly  from  Bellman’s  Principle  of  Optimality 
and  the  finiteness  of  the  Simplex  Method  for  linear  programming.  | 

BDP  requires  the  solution  of  a  great  number  of  linear  programs  in  Step 
1.  The  storage  requirements  may  be  prohibitive  for  a  large  or  detailed 
state  space,  Yt.  BDP  does,  however,  have  the  advantage  of  solving  small 
subproblems  and  of  following  a  single  pass  to  optimality  recursively  from 
period  t  back  to  1 .  The  method  may  even  be  implementable  when  one  has  a 
great  deal  of  advance  information  about  the  admissible  states  in  the  solution. 
For  more  general  problems,  however,  more  efficient  ways  to  characterize  the 
state  space  are  necessary.  In  the  next  section,  we  show  that  the  algorithms, 
PCSDLP  and  NDSDLP,  of  Chapters  III  and  IV  are  such  methods,  employing 
approximations  for  a  continuous  state  space. 

3.  Relation  to  the  Nested  Decomposition  Method 

In  this  analysis,  we  consider  the  optimization  problem,  DP(t),  at  some 
time  t : 


Zt(yt— i)  =  E&fcUfe— !>&)]'  (6) 

where  zt(yt-i,  6)  =  min  Vtzt((yt- i,  ft),  yt)  and 
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1,  St)  =  min  Ctzt  +  zt+  i(yt) 
subject  to 

Atzt=  yt- 1  +  St,  (7) 

—Btxt=  yt, 
xt>  0. 

Here,  we  would  like  to  have  an  explicit  representation  of  Zt+i{yt),  which 
we  approximated  with  the  quantizations  above.  One  way  to  find  zt+i{yt)  is 
to  exploit  its  convexity  properties.  The  nested  decomposition  and  piecewise 
approaches  do  this,  as  we  show  below. 

For  the  nested  decomposition  approach,  the  following  problem  at  stage 
t  is  solved  instead  of  DP(t)  : 


min  ctzt  -\-8t 
subject  to 

AtXt  =  yt— i  -f-  St, 

*«)*«+•«>  k  =  { ND(t)(a )) 

(— oli+iBt)xt  >  (alt+1at+i),l—  1, . . .  ,q,(ND{t){b)) 

xt>  0, 


{ND(t)) 


where  the  constraints  ND(t)(b)  keep  zt  feasible  and  ND(t)(a)  forms  an  outer 
linearization  of  the  convex  function,  Q(xt),  as  we  showed  in  Chapter  3.  The 
following  lemma  then  shows  the  equivalence  of  DP(t)  and  ND(t). 


Lemma  1.  For  all  t,  1  <  t  <  T,  x\  is  the  optimal  solution  of  zt{yt—  i,St)  if 
and  only  if  zj  solves: 

min  CfZtH-G(zt) 


subject  to  Atxt=  yt—  i  +  St, 


(8) 


Proof.  For  t  =  T,  we  have 


i)  =  E[Qr(*r— i,  £r)] 

=  E^t[  min  ctxt\Atxt  —  £r  4"  Bt—  i*r—  1, >  0]  (9) 

=  zr(Sr—i*r—i)- 

At  t  =  T  —  1, 


Qt—  i(*r— a,  £r— i)  ==  min  cr— i*r— i  4-  Qr(*r— i)  (10) 

subject  to 


At—iXt—i—  £t—i  +  By— a*r— a, 
*T— 1>  0, 


which,  by  (9),  is  equivalent  to  : 


min  ct—iZt—i  4*  Zt{Vt—i) 
subject  to 

At—  i*r— 1=  £r—  i  4*  Bjv_ a*r— a>  (H) 

Bt—iZt—i=  Vt—i, 

Zt — 1  ^  0. 

Therefore,  the  functions  zt—  i{vt— a>  (r—i)  and  Qt—i(zt—2>  (r—i)  are 
identical  and  have  correspondingly  the  same  optimal  solutions,  z*T_ v  The 
result  follows  by  induction  on  t.| 

Lemma  1  shows  that  the  cutting  plane  method  used  in  the  nested  decom¬ 
position  approach  is  equivalently  a  method  for  finding  the  function,  zt(yt— i), 
in  DP(t).  In  doing  this,  the  constraints  in  ND(t)(a)  linearise  the  convex  func¬ 
tion  in  the  neighborhood  of  yt—\. 

In  other  words,  at  stage  1 4- 1,  for  a  given  y\,  we  find 
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zt+i(vi)  —  E<tlzt+i(vii(t+i)] 

=  E(tl*t+i((t+i)(v}  +  (t+i)], 

=  *1+ 1 vt  (4+11 
~  (*i+iBt)xi  +  (4+i- 

By  convexity  and  for  (irJ+1Bt)it  +  a  support  of  zt+i(B*tt),  we 
obtain 


zt+i(BiZt)  >  (yr -f-  4+l-  (13) 

Now,  we  write  DP(t)  as 

min  ctxt  +  fit  ( DPP(t )) 

Atxt=  Vt-i  +  (t,(DPP(t)(a)) 

*t+i (Btxt)<  0t,  ( DPP(t)(b )) 

*t>  0, 

and  observe  that,  if  (zt,  fit)  is  such  that 

fit  <  (*i+iBt)zt  -f  p{+v  (14) 

then  (xt,  fit)  is  infeasible  in  DPP(t).  This  corresponds  (see  Figure  1.)  to 
finding  a  point  within  the  feasibility  space  of  the  previously  generated  cutting 
planes  that  is  not  feasible  for  Zt+t-  Another  plane  is  added  at  the  new  xt 
value  and  DPP(t)  is  solved  again.  The  process  is  repeated  until  (*tV?)  is 
found  such  that  0 J  =  2t+i(i?tzJ).  The  following  lemma  results: 

Lemma  2.  If  (jcJ,0J)  is  an  optimal  solution  to  ND(t)  and  0J  =  zt+\(Btx\), 
then  x\  is  an  optimal  solution  of  DP(t). 

Proof,  (x*,  0^)  is  optimal  for  the  reduced  constraints  in  ND(t)  and  x\  is 
feasible  for  DP(t)  since  fif  =  at+1(BtxJ),  therefore,  z\  is  optimal  for  DP(t).| 
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From  this  lemma,  NDSDLP  can  be  seen  as  a  method  for  an 

optimal  solution  to  the  dynamic  program,  DP(t).  We  further  observe  that 
the  discreteness  of  the  distribution  was  not  necessary  for  this  development. 
The  convexity  of  xt+i(yt)  is  sufficient  for  the  outer  linearization  to  properly 
bound  the  objective  function.  The  next  section  describes  how  PCSDLP  uses 
the  piecewise  linearity  of  zt+x(yt)  to  make  local  definitions  of  the  dynamic 
programming  recursion. 

4.  Relation  to  the  Piecewise  Method 

For  =  {£ J,  £2,  ,  },  from  Lemma  1,  it  follows  that 


k 

i=i 


t 


Figure  1.  Outer  linearizations. 


108 


Hence,  zt  end  Q  have  the  same  properties  as  functions  of  xt.  This 
implies  that  zt  is  piecewise  linear  in  xit  as  we  showed  for  Q(xt)  in  Chapter 
4.  Therefore,  we  can  write  DP(t)  as 

min  etxt  +  Sy— t)~ 

+jr,+a(B,+1[(^f')_,(fl.*.+f!+,)])  irvm 

subject  to 

Atxt=  yt—i  +  ft, 

=  1 . *, 

xt>  0. 

We  proceed  as  in  Chapter  4  to  write  PW(t)  as 

min  Ct^tH-*1 

subject  to 

Atit—  yt—i  +  ft,  {FWW[t)) 

—B[xt>  ft+i\;  =  l,...,*, 
xt>  0, 

where  ft*  is  a  constant,  and  ct,  B],  and  ft+i^  are  the  values  relative  to  the 
subproblem  bases  as  in  (3.12)  and  (3.15). 

Here,  the  function  zt((l/t—  i,  ft),!/t)  is  limited  to  the  linear  piece,  for 
which,  —B^xt  >  ft+iJ  for  all  j.  The  piecewise  algorithm  proceeds  by  look¬ 
ing  for  improvements  from  z*,  the  Ath  optimal  solution  found  to  PWW(t). 

A  sequence, 


*t((yt-i,  ft), v})  >  *t((yt- i,  ft), !/?)>•••>  **((Vt- 1, &)» vt),  (16) 

is  found,  such  that,  if  *t((yt—  i,  ft),!/?"*"1)  =  *t((Vt—  i,  ft),!/*)  f°r  feasible 
directions  from  yf,  then  the  x\  corresponding  to  B\x)  =  y)  is  optimal  for 
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DP(t),  and  we  can  proceed  to  find  *t(y<— i)-  (See  Figure  2,  and  note  that  the 
optimum  is  determined  for  zt  only,  so  that  Zt-f-i  need  not  be  at  a  minimum.) 

From  this  development,  we  hare 

Lemma  S.  If  the  piecewise  linear  approach,  PCSDLP,  terminates  with 

zt((yt—  i,  ft),  y*+1)  =  zt((yt_i,  ft),  yj),  then  the  associated  x*  is  an  optimal 
solution  to  zt((t/t—  1,  ft))  in  DP(t). 

4.  Conclusion 

In  Sections  3  and  4,  we  showed  that  the  piecewise  and  nested  decom¬ 
position  methods  found  values  of  xt  that  minimize  the  function  zt  for  a  given 
state  yt— i-  The  following  theorem,  which  follows  directly  from  Lemmas  3 
and  4,  states  this  result. 

Theorem  1.  The  nested  decomposition  and  piecewise  methods  presented  in 
Chapters  3  and  4  obtain  optimal  values  of  zt(yt— i,  ft)  at  each  stage  t  of  the 


Figure  2.  At  yf,  *t((yt-i,ft),y?)  = 
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dynamic  programming  procedure. 

PCSDLP  and  NDSDLP,  when  applied  to  this  backward  procedure,  must 
each  start  from  some  given  state,  yt.  The  dynamic  programming  approach 
of  enumerating  these  states  may  be  valuable  here.  If  we  let  yt  take  on 
representative  values  and  evaluate  each  of  the  associated  zt{yt—  1)  separately 
before  proceeding  to  stage  t  —  1,  some  computational  effort  may  be  saved  in 
both  the  nested  decomposition  and  piecewise  approaches. 


Ill 


CHAPTER  VH 


Computational  Results 
and 

Conclusions 


1.  Introduction 

In  this  chapter,  we  synthesize  our  previous  development  of  alternative 
methods  for  reducing  the  complexity  of  problems  under  uncertainty.  We  also 
present  some  results  from  solving  these  stochastic  dynamic  linear  programs. 
In  this  presentation,  we  wish  to  show  that  the  difficulties  in  solving  complex 
stochastic  programs  are  not  so  great  that  one  should  ignore  uncertainties. 
We  will  demonstrate  that  the  stochastic  solution  can  be  evaluated  without 
prohibitive  complications  and  that  our  techniques  may  prove  beneficial  in 
this  evaluation. 

Our  strategy  has  been  first  to  consider  deterministic  problems  and  then 
to  extend  them  to  the  SDLP  form.  Our  initial  step  in  this  approach  dealt  with 
the  properties  of  the  basis  for  different  scenarios.  In  Chapter  I,  we  showed 
that  this  analysis  may  result  in  finding  an  optimal  deterministic  solution. 
Beyond  this,  in  Chapter  2,  we  showed  that  a  bound  on  the  value  of  the 
stochastic  solution  could  be  obtained  before  proceeding  to  solve  SDLP.  We 
next  gave  three  algorithms  for  solving  this  problem.  In  the  following  sections, 
we  present  the  computational  properties  of  these  algorithms  in  solving  some 
examples  of  SDLP. 

Section  2  describes  our  computational  results  on  representative  test  prob¬ 
lems  and  compares  the  algorithms’  performance  with  standard  linear  pro- 
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gramming.  In  evaluating  these  methods,  we  also  consider  their  usefulness 
in  models  of  actual  phenomena.  In  Section  3,  we  discuss  areas  in  which 
stochastic  linear  programming  has  been  applied,  and  we  show  that  the  solu¬ 
tion  of  uncertainty  models  can  lead  to  prudent  decision-making.  Following 
this  analysis,  in  Section  4,  we  state  our  conclusions  on  solving  stochastic 
linear  programs  and  present  directions  for  future  research  in  this  area. 

2.  Computational  Results 

The  algorithms  of  Chapters  III,  IV,  and  V  have  been  programmed  in 
FORTRAN  on  the  SCORE  DECSystem  20  computer  at  Stanford  University. 
This  system  is  useful  for  observing  the  properties  of  algorithms,  but  it  is  not 
designed  for  solving  the  extremely  large  problems  modeled  in  some  stochastic 
dynamic  linear  programs.  Our  goal  was  not  to  draw  definitive  conclusions 
about  the  performance  of  each  algorithm,  but,  instead,  to  observe  some  of 
the  algorithms’  properties,  so  that,  on  systems  with  larger  processors,  they 
may  be  implemented  with  a  better  understanding  of  their  capabilities. 

Since  current  linear  programming  packages  take  great  advantage  of  the 
sparsity  (the  relative  number  of  zero  entries  in  the  coefficient  matrix),  we 
compared  our  algorithms’  performance  with  that  of  LPM-1,  a  program  writ¬ 
ten  by  J.  A.  Tomlin  and  revised  by  G.  Kochman  at  the  Systems  Optimization 
Laboratory.  LPM-1  employs  compact  storage  of  the  non-zero  entries  in  the 
coefficient  matrix,  performs  an  LU-decomposition  of  the  basis,  retains  the 
basis  inverse  in  product  form,  and  uses  a  merit  counting  sort  for  maintaining 
sparsity  in  the  inverse  (see  Pfefferkorn  and  Tomlin  [49]).  This  procedure  has 
proved  efficient  in  solving  linear  programs  because  the  coefficient  matrices  are 
characteristically  very  sparse.  Since  staircase  linear  programs  (and  SDLP’s) 
have  considerable  sparsity  in  their  structure,  it  has  been  conjectured  [47] 
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that  efficient  handling  of  these  elements  implicitly  uses  the  structure  and 
that,  except  on  very  large  problems  beyond  the  capacity  of  standard  linear 
programming  codes,  no  improvement  can  be  made  upon  the  packaged  codes. 
Our  procedures  have  an  advantage  over  these  codes  because  they  only  in¬ 
volve  small  sections  of  the  program  during  an  optimization.  They  do  not  re¬ 
quire  storage  of  the  entire  program  parameters  and,  thus,  may  handle  larger 
problems.  Apart  from  this  advantage,  we  want  to  observe  our  algorithms’ 
performance  on  smaller  problems  that  simplex  codes  handle  easily. 

The  nested  decomposition  method,  NDSDLP,  was  programmed  in  FOR¬ 
TRAN  as  NDST1.  This  program  includes  the  basic  algorithm  presented  in 
Chapter  III,  but  does  not  include  modifications,  such  as  column  passing.  For 
solving  individual  linear  programs  at  each  node,  NDST1  uses  the  procedures 
of  LPM-1.  It  saves  lists  of  the  current  basic  variables  at  each  node,  but  it 
does  not  maintain  the  LU-decomposition  for  each  node.  Instead,  it  keeps 
only  one  ETA  file  (the  elementary  matrices  in  the  product  form  of  the  basis 
inverse)  for  the  current  node  and  reinverts  the  basis  each  time  it  encounters  a 
new  node.  NDST1  also  includes  lower  and  upper  bound  capabilities  on  each 
variable  and,  hence,  does  not  require  Step  2'  of  NDSDLP. 

The  piecewise  approach,  PCSDLP,  was  programmed  as  PCST1.  It  in¬ 
cludes  the  feasibility  routines  of  NDSTl  for  finding  a  primal  feasible  solution 
of  SDLP,  and  it  follows  the  procedures  of  PCSDLP  given  in  Chapter  IV. 
rL  uxtu,  hxuwuvui ,  brave  vnyMtitiaat  -xavn/g  ’into  ifi  'xtfupnWnmi  Vaxex. 
Alternatively,  it  considers  each  subproblem  basis  individually  and,  therefore, 
may  be  forced  to  solve  larger  than  necessary  master  problems.  PCST1  also 
uses  lower  and  upper  bounds  on  all  variables  and  checks  for  possible  upper 
bound  violations  of  the  equation  (IV. 4a),  as  we  mentioned  above.  For  the 
case  of  artificial  variables  in  the  subproblem  basic  sets,  the  upper  bounds  are 
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always  included  as  additional  constraints. 


Storage  in  PCST1  is  larger  than  that  in  NDST1  because  of  the  need 
to  maintain  both  the  last  solution  for  and  the  solution  for  the  current 
auxiliary  problem  of  For  this  reason,  ETA  files  are  maintained  for  both 
problems. 


LBSMPX  was  implemented  in  the  program  LBS1.  This  program  uses 
the  basic  procedures  of  LPM-1  for  the  artificial  basis,  with  updated  transfor¬ 
mations  for  the  surplus  basis  blocks.  It  follows  the  cyclical  pricing  strategy 
described  in  Chapter  V  and  uses  a  “most  negative”  pricing  criterion  in  each 
period.  Pricing  is  restricted  to  that  period  as  long  as  the  least  reduced  cost  is 
negative.  Updates  of  the  surplus  basis  blocks,  Q\,  are  performed  according 
to  Ul,  U2,  U3,  and  U4  in  Chapter  V,  and  the  inverse  of  this  matrix  is  stored 
explicitly. 


We  have  chosen  a  set  of  representative  problems  to  show  the  algorithms’ 
computational  properties.  The  data  for  the  test  cases  appear  in  Table  1. 
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TABLE  1 


Program  Parameters 


Number 

of 

Conftraint* 

Number 

of 

Variable* 

Number 

of 

Node* 

Deniity 

(m) 

(n) 

00 

(/) 

PAl 

10 

22 

3 

.15 

RD1 

10 

25 

3 

.25 

NG1 

19 

35 

3 

.28 

PA2 

22 

50 

7 

.00 

RD2 

22 

57 

7 

.14 

PA3 

46 

106 

15 

.04 

RD3 

46 

121 

15 

.07 

*  (number  of  non-iero  element!) /(total  number  of  element*) 


These  examples  -were  chosen  to  represent  the  types  of  problems  found  in 
applications. 

The  first  cases,  RD1,  RD2,  and  RD3,  were  selected  from  a  control 
problem  in  Davis  [20]  and  modified  to  reflect  a  production  planning  model 
as  in  Beale  et  a/  [7].  In  this  format,  the  expected  present  value  of  profits  is 
maximized  for  a  firm  planning  the  production  and  storage  of  a  number  of 
products  that  require  common  resources.  Their  sales  are  determined  by  an 
uncertain  demand.  The  decision  variables  in  this  model  are  : 

z»,t  —  the  amount  of  product  i  sold  in  period  t, 

yi't  —  the  amount  of  i  produced  in  period  t, 

8^  —  the  amount  of  product  t  in  stock  at  the  end  of  period  t. 

The  linear  program  for  this  problem  is 
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(1) 


max  Et^l  J2i=l(Pi,txi,t  ki,t*i,t) 

8.  t.  32lmi  Vi,t<  Rt>  for  aW  *> 

«i,t-i  +  Vi,t  —  2*,t=  Si't,  for  all  t  and  *, 

0  <  xi,t  <  d»,t,  for  all  £  and  t, 

where  Rt  is  the  total  capacity  in  period  £,  ptjt  is  the  selling  price  of  product 
i,  Ci,t  is  the  production  cost  of  product  t,  k^t  is  the  holding  cost  of  stock  in  t, 
T  is  the  planning  horizon,  and  q  is  the  number  of  products.  In  this  example, 
we  included  random  uncertainties  in  di,t  and  formulated  (1)  as  in  SDLP  for 
optimizing  the  expected  value  of  profits  over  the  planning  horizon. 

The  examples,  PA1,  PA2,  and  PA3,  were  contributed  by  P.  Abrahamson 
[1]  and  represent  the  discrete  approximation  of  a  continuous  time  linear 
program.  This  program  represents  the  optimal  control  of  a  particle  subject 
to  state  space  constraints.  It  is  written  as 

min  /o  u(t)eatdt 

s.t.  x(t)  —  —  z(s))ds=  1, 

“(f)  <  1,  (2) 

*(f)<  .7  — .2 £, 

z(f),“(f)>  0, 

where  we  considered  uncertainty  in  the  state  to  state  transitions  for  x[t). 
The  application  of  uncertainty  to  continuous  problems,  as  suggested  by  this 
model,  is  an  important  area  for  possible  future  applications. 

The  last  example,  NG1,  is  derived  from  the  exhaustible  resource  model  of 
Chapter  I.  It  includes  uncertainty  in  investment  in  different  technologies.  This 
example  is  included,  because  the  basis  must  be  changed  significantly  from  the 
myopic  first  period  solution  to  the  optimal  stochastic  solution.  This  process 
requires  more  surplus  basic  variables  in  the  first  period  and  demonstrates 
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TABLE  2 


Sample  Timet  and  Iteration » 


LPM-1  NDST1 _ PCST1 _ UBS1 


Time 

(CPU.) 

No. 

of 

Iter. 

Time 

(CPU.) 

No. 

of 

Iter. 

Time 

(CPU.) 

No. 

of 

Iter. 

Time 

(CPU*) 

No. 

of 

Iter. 

PA1 

.18 

12 

.14 

12 

.10 

12 

.20 

11 

RD1 

.30 

15 

.23 

12 

.42 

13 

.33 

15 

NG1 

.48 

18 

.83 

20 

1.01 

23 

.56 

17 

PA2 

.04 

20 

1.71 

20 

2.41 

30 

1.05 

28 

RD2 

1.50 

31 

.01 

32 

1.34 

33 

1.73 

31 

PAS 

4.01 

61 

3.30 

88 

5.61 

70 

4.06 

64 

RD3 

7.13 

85 

4.80 

67 

7.64 

71 

6.72 

50 

some  of  the  complications  of  stochastic  programs. 

The  solution  times  and  number  of  simplex  iterations  required  for  the 
optimal  solution  of  these  problems  are  presented  in  Table  2.  The  times 
are  in  CPU  seconds  for  solutions  from  a  “cold  start”  (an  infeasible  basis), 
and  represent  computing  time  excluding  time  for  data  input.  The  results 
show  that,  on  small  programs,  our  algorithms  may  be  competitive  with 
codes  using  sparsity  techniques  alone,  such  as  LPM-1.  This  is  especially 
significant  because  NDSTl,  PCST1,  and  LBS1  are  experimental  codes  and 
are  not  programmed  for  maximum  computing  efficiency.  Small  problems, 
however,  are  not  representative  of  actual  applications,  and  experience  with 
larger  programs  is  mandatory  for  true  comparisons.  Again,  we  do  not  intend 
to  prove  the  superiority  of  our  methods,  but  only  to  show  how  they  perform. 
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TABLE  3 


Sample  Times  per  Iteration 

LPM-1 
Time/No.  of 
Iteration* 

NDST1 
Time/No.  of 
Iteration* 

PCST1 
Time/No.  of 
Iteration* 

LBS1 

Time/No.  of 
Iteration* 

PAX 

.015 

.012 

.016 

.018 

RD1 

.020 

.018 

.032 

.022 

NG1 

.027 

.028 

.044 

.033 

PA2 

.032 

.058 

.080 

.037 

RD2 

.051 

.028 

.041 

.056 

PA3 

.066 

.048 

.071 

.078 

RD3 

.110 

.073 

.108 

.114 

NDST1  performs  best  overall  on  these  examples,  but,  where  it  performs 
poorly  we  can  observe  some  of  its  properties.  In  NG1,  a  large  number  of 
cuts  were  required  on  the  first  period  and  many  suboptimisations  were  per¬ 
formed,  slowing  its  convergence  to  an  optimum.  We  also  see  here  that  these 
suboptimizations  required  NDST1  to  perform  25  percent  more  iterations  than 
PCST1,  which  maintains  subproblem  optimality.  This  difficulty  also  appears 
in  PA2  where  the  added  effort  of  reinversion  at  each  pass  from  node  to  node 
led  to  NDSTl’s  having  a  larger  average  time  per  iteration  than  LPM-1  (see 
Table  3.).  We  note  that  the  small  growth  of  non-zeroes  in  the  basis  inverse 
of  PA2  is  a  good  example  of  LPM-l’s  advantage  in  performing  backward 
and  forward  transformations  very  quickly.  On  the  other  examples,  because 
NDST1  solves  the  smallest  problems,  it  requires  the  least  amount  of  time  per 
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iteration  despite  reinyersion. 

An  interesting  property  observed  in  the  solutions  by  NDST1  concerns  the 
set  of  basic  variables  in  the  master  problem  as  cutting  planes  are  added  from 
the  subproblems.  In  our  examples,  we  often  saw  that  one  set  of  basic  variables 
would  remain  constant  in  the  master  problem  from  iteration  to  iteration  and 
that  the  additional  basic  elements,  would  be  chosen  from  another  small  set 
of  variables.  The  algorithm  would  choose  a  member  of  the  additional  set  on 
one  iteration,  replace  it  with  another  member  in  the  next  iteration,  and  then 
bring  the  first  element  back  into  the  basis  at  master-subproblem  optimality. 
For  Xg,  the  constant  set  of  basic  variables,  and  {xi,  z2},  the  additional  basic 
variables,  the  basis  in  the  master  problem  would  follow  the  pattern: 

Iteration  Basic  Set 

1  {XB} 

2  {XB,zx} 

3  {xBlx2y 

4  { XB ,  *1,  x2}~  optimal. 

It  has  been  observed  [1]  that  the  deterministic  problem  often  brings 
additional  variables  into  the  basis  and  then  adjusts  them  in  subsequent 
iterations  without  dropping  them  from  the  basis.  This  observation  led  to 
the  column  passing  technique  we  discussed  above.  Our  observations  indicate 
that  the  set  of  columns  passed  forward  to  the  subproblems  should  include 
columns  that  were  in  the  basis  but  were  deleted.  This  would  allow  for  the 
entire  set  of  columns  to  obtain  the  optimal  weights  in  a  single  pass. 

The  behavior  we  have  observed  may  be  indicative  of  the  hedging  effect  in 
stochastic  problems.  In  our  example  above,  z\  may  be  included  to  optimise 
the  first  descendant  scenario  and  z2  may  correspond  to  optimising  with 


TABLE  4 


Number  of  Firet  Ptriod  Optimisation! 


NDST1 

PCSTl 

Number 

of 

P&iiee 

Number 

of 

Piuim 

PA1 

6 

4 

RD1 

8 

T 

NG1 

14 

11 

PA2 

8 

6 

RD2 

e 

7 

PA3 

12 

10 

RD3 

15 

13 

the  second  subproblem.  The  algorithm  tries  each  of  these  activities  before 
deciding  on  the  optimal  combination  in  iteration  4.  By  including  each  of  these 
potentially  banc  variables  in  the  set  corresponding  to  the  columns  passed 
forward,  we  may  avoid  eneessive  optimisation. 

PCST1  performed  consistently  worse  than  NDST1  based  on  CPU  time, 
but  its  strength  of  not  requiring  suboptimization  appears  in  Table  4.  This 
chart  displays  the  number  of  times  the  master  problem  in  period  1  was 
solved.  PCST1  requires  fewer  of  these  optimizations  and  passes  to  the 
subproblems.  It  thus  decreases  the  possibility  of  lengthy  suboptimization. 
PCST1  solves,  however,  a  larger  master  problem  than  NDST1,  because  it 

includes  the  degeneracies  associated  with  repeated  subproblem  bases.  The 
basis  list  approach  presented  in  Chapter  4  should  alleviate  some  of  the 
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difficulty  here. 

LBS1  also  maintains  more  information  than  is  necessary  and  its  use  of 
the  triangular  basis  differs  only  slightly  from  the  factorization  for  sparsity  in 
LPM-1.  LBSl’s  additional  bookkeeping  is  most  evident  on  smaller  programs 
and  begins  to  be  efficient  for  larger  problems  as  we  start  to  see  in  RD3. 

These  examples  presented  some  of  the  properties  of  NDST1,  PCST1,  and 
LBSl.  Knowledge  of  these  attributes  should  be  helpful  in  determining  what 
method  to  use  for  a  specific  problem.  To  determine  this,  as  we  stated  above, 
the  deterministic  scenarios  should  be  evaluated  first  and  then  the  stochastic 
problem  should  be  solved.  Because  of  its  simplicity,  we  would  recommend 
NDST1  as  a  first  method  to  try,  followed  by  PCST1,  if  convergence  to 
optimality  is  slowed  by  excessive  suboptimization.  LBSl’s  implementation  is 
more  determined  by  efficient  coding  and  its  evaluation  should  be  made  only 
once  its  procedures  have  been  efficiently  coded  for  larger  problems  on  faster 
processors. 

3.  Areas  of  Application 

We  have  mentioned  that  stochastic  dynamic  linear  programs  may  be 
found  in  many  real-world  contexts.  In  this  section,  we  discuss  some  actual 
applications  of  these  problems  in  decision-making.  We  mention  these  ex¬ 
amples  to  demonstrate  the  range  of  possible  implementations  of  our  tech¬ 
niques. 

The  first  use  of  stochastic  linear  programs  was  by  Dantzig  and  Ferguson 
[21]  on  a  problem  of  airline  scheduling  with  uncertain  passenger  demand. 
Using  a  modification  of  the  transportation  simplex  method,  they  demonstrated 
that  the  net  expected  costs  for  meeting  the  carrier’s  demand  could  be  reduced 
from  $1,666,000  to  $1,524,000  by  considering  the  complete  distribution  of 
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demand  instead  of  only  the  expected  values  of  those  demands. 

More  recent  uses  of  stochastic  linear  programs  include  Manne’s  analysis 
in  [44]  of  the  decision  to  develop  breeder  fission  technology.  Another  analysis 
of  nuclear  fuel  choices  appears  in  Avi-Itzhak  and  Connally  [4],  wherein, 
they  use  the  scenario  approach  of  assuming  decisions  based  on  deterministic 
solutions  and  evaluating  the  expected  cost  of  these  decisions  under  different 
future  scenarios. 

Other  areas  of  application  are  found  in  the  investment  management  of 
bank  portfolios.  Aghili  et  ai  [2]  evaluated  the  decisions  of  a  small  midwestern 
bank  in  allocating  their  assets  among  a  set  of  securities  and  loan  positions. 
They  used  the  bank  management’s  criteria  for  constraints  and  formulated  a 
two-stage  stochastic  linear  program  with  uncertain  interest  rates  and  bond 
yields.  Their  results  using  only  financial  and  accounting  constraints  reflected 
the  extreme  point  property  of  optimal  linear  program  solutions.  Their  model 
would  have  placed  almost  all  of  the  bank’s  portfolio  into  municipal  bonds, 
forgoing  all  commercial  loans  and  mortgages. 

The  difficulty  Aghili  et  a 1,  encountered  arose  from  their  use  of  the  risk 
neutral  linear  objective  function.  To  improve  their  results,  they  incorporated 
management’s  concerns  into  additional  policy  constraints  and  obtained  solu¬ 
tions  that  were  much  closer  to  the  bank’s  own  decisions.  From  this  analysis, 
management  became  aware  of  the  need  to  consider  different  future  economic 
environments.  Previously,  they  had  used  only  one  forecast  for  the  future 
but  were  now  able  to  consider  several  scenarios.  This  openness  to  changing 

money  conditions  gave  me  Danx  management  more  options  to  consiaer  in 
their  policy  decisions  and  led  them  to  reshape  some  of  their  ideas  on  the 
bank’s  operations. 

The  bank  balance  sheet  model  is  a  good  example  of  decision  makers’ 


use  of  a  stochastic  model.  Difficulties  often  arise  in  the  implementation  of  a 
model’s  optimal  solution,  but,  in  Aghili  et  a/’s  experience,  the  model  served  to 
inform  management,  and  it  became  an  important  policy  instrument.  Another 
example  of  this  utility  of  stochastic  linear  models  appears  in  Gaither’s  study 
[24]  of  commercial  fishing  seasons.  That  research  demonstrated  that  sub¬ 
stantial  savings  could  be  accrued  by  proper  fishing  regulation  in  Alaska,  but 
modifying  those  regulations  required  difficult  policy  decisions  that  would  not 
be  immediately  forthcoming.  Despite  these  problems,  Gaither’s  model  also 
played  an  important  role  in  demonstrating  to  the  decision  makers  how  they 
should  structure  their  strategy  and  what  concerns  they  should  have. 

4.  Conclusion 

We  have  presented  several  approaches  for  solving  linear  models  that  in¬ 
volve  decisions  made  over  time  in  uncertain  environments.  We  have  showed 
that  the  resulting  stochastic  dynamic  linear  programs  can  be  analyzed  through 
the  use  of  deterministic  scenario  solutions  and  the  subproblem  solution  methods 
of  nested  decomposition,  piecewise  path  following,  and  local  basis  factoriza¬ 
tion.  This  analysis  rests  upon  the  fundamental  properties  of  the  basis  in 
linear  programming,  the  understanding  of  which  is  crucial  in  understanding 
the  characteristics  of  a  stochastic  solution. 

The  stochastic  solution  succeeds  because  it  allows  for  a  regularization 
of  the  extreme  point  solutions  that  occur  in  different  deterministic  scenarios. 
This  property  should  be  carefully  considered  in  evaluating  any  stochastic 
solution  because  it  can  be  a  powerful  stabilizing  force  in  determining  optimal 
decisions.  Within  the  framework  of  the  stochastic  model,  a  single  decision 
can  be  chosen  that  will  reflect  each  of  the  future  scenarios  included  in  the 


The  methods  we  presented  for  solving  these  SDLP’s  demonstrated  credible 
performance  on  smaller  examples.  These  results  are  encouraging  for  future 
implementations  on  large-scale  problems.  This  experimentation  on  large 
models  should  be  the  major  direction  for  future  research. 

The  true  test  of  our  algorithms’  benefits  and  of  the  usefulness  of  the 
SDLP  program,  in  general,  must  come  from  actual  policy  models.  Our 
presentation  of  examples  from  production  planning,  exhaustible  resources, 
and  portfolio  management  demonstrate  that  many  applications  of  SDLP  are 
possible.  The  inclusion  of  uncertainty  into  a  linear  model  extends  the  range 
of  decision  factors  and  helps  stabilize  policy  strategies  in  dynamic  settings. 
These  attributes  may  bring  significant  gains  to  the  use  of  models  in  decision¬ 
making. 

Our  primary  goal  has  been  to  show  that,  by  exploiting  the  stochastic 
linear  program’s  special  structure,  models  may  be  able  to  include  uncertain 
parameters  without  becoming  hopelessly  complex.  By  proceeding  from  deter¬ 
ministic  scenarios  through  the  stochastic  solution,  the  modeler  can  be  aware 
of  the  value  of  proceeding  to  more  complex  stages  and  can  develop  a  better 
understanding  of  the  nature  of  the  deterministic  linear  programming  solu¬ 
tion  and  its  relation  to  the  stochastic  strategy.  To  complete  this  analysis 
most  effectively,  the  computational  costs  of  the  stochastic  program  must  be 
reduced.  We  hope  that  our  work  will  help  make  these  reductions  possible, 
thereby,  enabling  models  to  become  more  powerful  tools  in  affecting  current 
decisions  in  an  uncertain  world. 
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SOL  80-29:  Solution  Methods  for  Stochastic  Dynamic  Linear  Programs, 
f  John  R.  Blrge 

l 

'  Linear  programs  have  been  formulated  for  many  practical  situations 
that  require  decisions  made  periodically  through  time.  These  dynamic 
linear  programs  often  involve  uncertainties .  Deterministic  solutions 
of  these  problems  may  lead  to  costly  incorrect  decisions,  and,  when  a 
stochastic  solution  is  attempted  the  problem  may  become  too  large.  In 
this  report,  we  present  methods  for  reducing  the  computational  cost  of 
these  stochastic  programs,  and  we  show  conditions  under  which  the 
stochastic  program  need  not  be  solved.  Our  methods  are  based  on  the 
large-scale  programming  techniques  of  decomposition,  partitioning,  and 
basis  factorization.,. 
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